library(ggplot2)
library(pander)
library(knitr)
library(tidyr)
library(mgcv)
library(diversitree)
library(repmis)
library(lme4)
library(ape)
library(geiger)
library(grid)
library(gridExtra)
library(nlme)
library(phytools)
library(brms)
library(ggridges)
library(caper)
library(purrr)
library(reshape2)
library(ggExtra)
library(grid)
library(dplyr)
library(BAMMtools)
library(gridExtra)
library(car)
library(coda)
library(MuMIn)
library(parallel)
library(HDInterval)
library(phangorn)
library(kableExtra)
library(RColorBrewer)
library(tibble)
library(stargazer)
library(ggraph)
#devtools::install_github("Ax3man/phylopath")
library(phylopath)
library("EBImage") # for images
library("ggtree")
source("bamm_extinction/functions/check_and_fix_ultrametric.R")
source("functions/essim.R")
source("functions/diversification_rate_calculation.R")
source("functions/Vdodge_function.R")
source('functions/symmetrise_scale.R')
source('functions/Correlation_matrix.R')
source('functions/gheatmap2.R') #Gets rid of spaces between tiles

Here we provide the supplementary material for our study on sexual selection and extinction risk in passerine birds.

Compiling the data

For this project we aimed to compile data on a large set of species in the Passerine bird order that would allow us to account for many evolutionary and environmental predictors of speciation. The following details the source and use of the data.

Speciation Rate

To obtain speciation rate estimates we used Bayesian Analysis of Macroevolutionary Mixtures (BAMM). Additionally, we compared our BAMM output files (event data) against a pre-existing extinction rate dataset (Harvey et al. 2017).

Rather than obtaining one estimate of speciation rate from one tree, we ran BAMM 100 times on 100 trees plus an MCC tree to obtain uncertainty estimates of the speciation rate generated from the phylogeny created by Jetz et al. (2012). We also estimated extinction rate from BAMM.

Sexual Selection

Our proxy for sexual selection was based on measures of sexual dichromatism. Notably, we use two existing datasets; a complete set of male and female plumage scores in Passerines from Dale et al. (2015) and reflectance data from 1000 birds from Armenta et al. (2008).

Briefly, Dale et al. (2015) obtained male and female plumage score, which can be divided to obtain a sexual dichromatism index (SDi). The plumage score was scored using the Handbook of the Birds of the World Del Hoyo et al. (2011). By scanning images of males and females in this book across multiple patches Dale et al. (2015) were able to obtain mean plumage scores from RGB values. Notably concern about the accuracy of these measurements were raised, thus this dataset was previously compared with a more thorough reflectance dataset of Australian birds, the two datasets correlated but was expectedly noisy. Here we compare the Dale et al. (2015) data against another reflectance dataset from Armenta et al. (2008) using colour discriminability as the parameter.

Alongside sexual dichromatism data there is also partial data for other metrics analysed in Dale et al. (2015) such as: Body size, tropical life history, Sexual selection, Cooperative breeding and Migration. We use a phylogenetic PCA (PC1) that estimates male-bias sexual selection from a species level of sexual size dimorphism, social polygyny and paternal care.

Environmental predictors

We obtained species range distribution maps from Birdlife International (through IUCN) (BirdLife International and Handbook of the Birds of the World 2017). These species range maps cover nearly all species of birds. From these spatial information files we were able to obtain estimates of the following (note that not all the data extracted was subsequently used in the analysis):

  • Range Size
  • Average and standard deviation in 19 bioclimatic variables (each range randomly sampled 1000 times)
  • Average and standard deviation in 19 bioclimatic variables from the last-glacial maximum (LGM) and the last-inter-glacial (LIG)
  • Net primary productivity (NPP) estimates and variability
  • Average and standard deviation of human population density in the species ranges.

The code for extracting the environmental data from the species is provided in another R-markdown document Github link, however no raw data is provided alongside this file due to the large file size of the shapefile and raster information. We have put the extracted environmental variables in a csv filed called complete.dataframe.csv for convenience. Briefly for each of the ~ 6,000 species, we:

1) Randomly sampled the range polygon 1000 times:

#Increse the iterations (defaukt is 4) so we can obtain complete samples of each range
bird.points <- lapply(bird.ranges, 
                       function(x) {spsample(x, n=1000, type="random", iter = 30)})

2) Extracted the bioclim or other (e.g. NPP) data from each of the points:

bird.values <- lapply(bird.points, function(x) {raster::extract(bioclim_data, x)})
bird.values <- as.list(data.frame((bird.values)))

3) Summarised the extracted data across the 1,000 points into a summary value of interest (means and se) for each species and exported that data.

#Obtain means
bird.frame <- lapply(bird.values, function(x) {as.data.frame(x)})
bird.summary <- lapply(bird.frame, function(x) {
(as.data.frame(apply(x,2,mean, na.rm =T )))})

#transpose
bird.means <- t(as.data.frame(bird.summary))
bird.means <- (split(bird.means, 1:19))

#Now add column of species name: Same order carried through
bird.means <- cbind.data.frame(shps.jetz, bird.means)

rownames(bird.means) <- NULL
colnames(bird.means) <- c("binomial","bioclim1", "bioclim2", "bioclim3", "bioclim4", "bioclim5", "bioclim6", "bioclim7", "bioclim8", "bioclim9", "bioclim10", "bioclim11", "bioclim12", "bioclim13", "bioclim14", "bioclim15", "bioclim16", "bioclim17", "bioclim18", "bioclim19")

#Write csv
write.csv(bird.means, 'data/bird.means.csv')

4) This process was repeated when we used gridded data for the LGM, LIG, NPP and human population density, while range size was obtained from the following code:

bird.range.size <- sapply(bird.ranges,
                       function(x) {(area(x))})
bird.range.size <- sapply(bird.range.size, function(x) {sum(x)})
bird.range.size <- as.data.frame(bird.range.size)
bird.range.size <- cbind.data.frame(shps.jetz, bird.range.size)
rownames(bird.range.size) <- NULL
colnames(bird.range.size) <- c("binomial", "range.size.m2")

#write.csv
write.csv(bird.range.size, 'data/bird.range.size.csv')

Generating biologically relevent predictors

Here we assess the relationship between sexual selection and extinction risk. However, in doing so we attempt to take into account as many other predictors of extinction as possible, primarily through environmental variables. Basically, our model structure seeks to contain:

Extinction/Diversification ~ Sexual selection
+ Range size
+ Short temporal variability of temperature (mean BIOCLIM4)
+ Spatial variability of temperature (PCA1) [residual.PCA1]
+ Long-term variability of temperature (LIG)
+ NPP

These predictors can be obtained from the compiled datasets seen here:

plumage.scores <- read.csv('data/plumage_scores.csv')
#Generate sexual dichromatism score: 
plumage.scores$SDi <- abs(plumage.scores$Male_plumage_score - plumage.scores$Female_plumage_score)

#Make DF with enviro variables and plumage scores
complete.dataframe <- read.csv('data/complete.dataframe.csv')
complete.dataframe <- left_join(plumage.scores %>% dplyr::select(binomial, Male_plumage_score, Female_plumage_score, SDi), complete.dataframe, by = "binomial")

Spatial Variability PCAs

To obtain estimates of spatial variability we can obtain PCAs of the standard errors of the bioclimatic variables (excluding seasonality in temperature and precipitation). From there we can run GAM models with the PCA components and range size and extract residuals. We do this as we expect variability will increase alongside increased range size, thus taking the residuals allows us to correct for this association.

Table S1: Loadings for the first two PCs, when we conduct a PCA on the variation (standard error, where n ~ 1,000) pf each bioclim variable except the standard error of temperature seasonality (se.bioclim4) and the standard error of precipitation seasonality (se.bioclim15), as these are already measures of variability and based on other bioclim variables. We use the standard error as, while we attempted to sample each species range 1,000 times, some points returned NA values and would thus affect the standard deviation, this inconvenience should not affect standard error.

restricted.data <- complete.dataframe %>% drop_na(bioclim1) #Drop species without environmental data
PCA.bioclim <- prcomp(restricted.data[c('se.bioclim1', 'se.bioclim2', 'se.bioclim3', 'se.bioclim5', 'se.bioclim6', 'se.bioclim7', 'se.bioclim8', 'se.bioclim9', 'se.bioclim10', 'se.bioclim11', 'se.bioclim12', 'se.bioclim13', 'se.bioclim14', 'se.bioclim16', 'se.bioclim17', 'se.bioclim18', 'se.bioclim19')], #Select se.bioclim
                      scale = TRUE, center = TRUE)
PCA.predictions <- predict(PCA.bioclim)
restricted.data <- cbind(restricted.data, PCA.predictions)

#Create table with PC1 and PC2
as.data.frame(PCA.bioclim$rotation[,1:2]) %>% `colnames<-`(c("PC1", "PC2")) %>% pander()
  PC1 PC2
se.bioclim1 -0.3328 0.08264
se.bioclim2 -0.2383 0.05458
se.bioclim3 -0.2317 -0.06584
se.bioclim5 -0.312 0.08091
se.bioclim6 -0.3288 0.08275
se.bioclim7 -0.2934 0.09748
se.bioclim8 -0.2916 0.1181
se.bioclim9 -0.3187 0.1125
se.bioclim10 -0.3111 0.06704
se.bioclim11 -0.3296 0.08775
se.bioclim12 -0.1386 -0.4246
se.bioclim13 -0.1596 -0.3083
se.bioclim14 -0.01133 -0.3857
se.bioclim16 -0.1708 -0.3254
se.bioclim17 -0.01923 -0.3932
se.bioclim18 -0.1433 -0.3066
se.bioclim19 -0.01797 -0.3815
#run GAM
PC1.gam <- mgcv::gam(PC1*(-1) ~ s(log(range.size.m2)), data = restricted.data, family = "gaussian") #Inverse as the PC1 loads to the negative (COUNTER-NTUITIVE)

#take residuals
restricted.data$residuals.PC1 <- residuals.gam(PC1.gam) 

#We can do the same for PC2
PC2.gam <- mgcv::gam(PC2 ~ s(log(range.size.m2)), data = restricted.data, family = "gaussian")

#Take residuals
restricted.data$residuals.PC2 <- residuals.gam(PC2.gam)

#plot
par(mfrow = c(1,2))
plot.gam(PC1.gam, residuals = T, main = 'PC1 residuls (temp)')
plot.gam(PC2.gam, residuals = T, main = 'PC2 residuals (precip)')

Figure S1: The relationship between spatial variability in the temperature components (PC1) and log-range size is relatively strong. But not as strong for spatial variability in precipitation (PC2). In the analysis we only used the residuals from PC1

Long term climate variability (LIG anomaly)

To gain estimates for change in climate over the past ~130,000 years we can use the difference in bioclim variables between the LIG and present values. The plots show the two PCs, broadly representing temperature and prescipitation. Data has also been provided for the last-glacial maximum (LGM) but was not used in the analysis.

Table S2: Loadings for the first two PCs of each of the PCAs for the difference in bioclimatic variables between today and the LIG. Here, PC1 is more heavily loaded for absoluted temperature difference, while PC2 is more heavily loaded for absolute difference in precipitation.

historical.variation.data <- as.data.frame(restricted.data[1])

#FOR LIG

historical.variation.data$bio1.LIG.diff <- abs(restricted.data$bioclim1 - restricted.data$LIG.bi1)
historical.variation.data$bio2.LIG.diff <- abs(restricted.data$bioclim2 - restricted.data$LIG.bi2)
historical.variation.data$bio3.LIG.diff <- abs(restricted.data$bioclim3 - restricted.data$LIG.bi3)
historical.variation.data$bio4.LIG.diff <- abs(restricted.data$bioclim4 - restricted.data$LIG.bi4)
historical.variation.data$bio5.LIG.diff <- abs(restricted.data$bioclim5 - restricted.data$LIG.bi5)
historical.variation.data$bio6.LIG.diff <- abs(restricted.data$bioclim6 - restricted.data$LIG.bi6)
historical.variation.data$bio7.LIG.diff <- abs(restricted.data$bioclim7 - restricted.data$LIG.bi7)
historical.variation.data$bio8.LIG.diff <- abs(restricted.data$bioclim8 - restricted.data$LIG.bi8)
historical.variation.data$bio9.LIG.diff <- abs(restricted.data$bioclim9 - restricted.data$LIG.bi9)
historical.variation.data$bio10.LIG.diff <- abs(restricted.data$bioclim10 - restricted.data$LIG.bi10)
historical.variation.data$bio11.LIG.diff <- abs(restricted.data$bioclim11 - restricted.data$LIG.bi11)
historical.variation.data$bio12.LIG.diff <- abs(restricted.data$bioclim12 - restricted.data$LIG.bi12)
historical.variation.data$bio13.LIG.diff <- abs(restricted.data$bioclim13 - restricted.data$LIG.bi13)
historical.variation.data$bio14.LIG.diff <- abs(restricted.data$bioclim14 - restricted.data$LIG.bi14)
historical.variation.data$bio15.LIG.diff <- abs(restricted.data$bioclim15 - restricted.data$LIG.bi15)
historical.variation.data$bio16.LIG.diff <- abs(restricted.data$bioclim16 - restricted.data$LIG.bi16)
historical.variation.data$bio17.LIG.diff <- abs(restricted.data$bioclim17 - restricted.data$LIG.bi17)
historical.variation.data$bio18.LIG.diff <- abs(restricted.data$bioclim18 - restricted.data$LIG.bi18)
historical.variation.data$bio19.LIG.diff <- abs(restricted.data$bioclim19 - restricted.data$LIG.bi19)

historical.variation.data <- historical.variation.data %>% drop_na(bio1.LIG.diff)

#Run a PCA of the difference removing the variation variables (4 and 15)

#For LIG
PCA.LIG.bioclim <- prcomp(historical.variation.data[c(
'bio1.LIG.diff',
'bio2.LIG.diff',
'bio3.LIG.diff',
'bio5.LIG.diff',
'bio6.LIG.diff',
'bio7.LIG.diff',
'bio8.LIG.diff',
'bio9.LIG.diff',
'bio10.LIG.diff',
'bio11.LIG.diff',
'bio12.LIG.diff',
'bio13.LIG.diff',
'bio14.LIG.diff',
'bio16.LIG.diff',
'bio17.LIG.diff',
'bio18.LIG.diff',
'bio19.LIG.diff')],
                      scale = TRUE, center = TRUE)

#Create table with PC1 and PC2
as.data.frame(PCA.LIG.bioclim$rotation[,1:2]) %>% `colnames<-`(c("PC1.LIG", "PC2.LIG")) %>% pander()
  PC1.LIG PC2.LIG
bio1.LIG.diff -0.2956 -0.01943
bio2.LIG.diff -0.2442 0.0745
bio3.LIG.diff -0.1436 0.2085
bio5.LIG.diff -0.2867 -0.2566
bio6.LIG.diff -0.4094 0.03346
bio7.LIG.diff -0.4012 -0.09706
bio8.LIG.diff 0.07108 -0.2385
bio9.LIG.diff -0.3365 -0.00753
bio10.LIG.diff -0.07014 -0.3766
bio11.LIG.diff -0.4201 -0.004
bio12.LIG.diff -0.05181 0.3337
bio13.LIG.diff -0.2552 0.2591
bio14.LIG.diff 0.05701 0.3298
bio16.LIG.diff -0.1946 0.3164
bio17.LIG.diff 0.03896 0.3603
bio18.LIG.diff 0.06173 0.277
bio19.LIG.diff 0.08743 0.2856
#Now we can predict the PCA results: 
PCA.LIG.predictions <- as.data.frame(predict(PCA.LIG.bioclim)*(-1)) #So that higher numbers mean more variation
PCA.LIG.predictions <- rename(PCA.LIG.predictions, PC1.LIG = PC1,
       PC2.LIG = PC2)

#Bind them to dataframe, taking only the first two PCAs
historical.variation.data <- cbind(historical.variation.data, PCA.LIG.predictions[1:2])

#Now to the restricted dataframe
restricted.data <- right_join(restricted.data, historical.variation.data %>% dplyr::select(binomial, PC1.LIG, PC2.LIG), by = 'binomial')

PCA.plot <- historical.variation.data %>% ggplot(aes(x = PC1.LIG, y= PC2.LIG))+
  geom_point(shape = 21)+
  theme_minimal()

PCA.plot.m <- ggExtra::ggMarginal(
  p = PCA.plot,
  type = 'density',
  margins = 'both',
  size = 5,
  colour = 'black',
  fill = 'gray'
)

grid.arrange(PCA.plot.m)

Figure S2: The distribution of the two PCA compontents for the absolute difference in bioclimatic variables between today and the LIG.

Correlations between environmental predictors

We can check the correlations between environmental predictors that we will use in our model. Specifically we want to know whether the following environmental variables are correlated:

  • Range size
  • Short temporal variability of temperature (mean BIOCLIM4)
  • Spatial variability of temperature (PCA1) [residual.PCA1]
  • Long-term variability of temperature (LIG)
  • NPP

To check some of the correlations between the environmental predictors to be used in the PGLS models we can inspect the following correlation plot with the correlation value plotted:

#Draw plot
restricted.data$log.range.size <- log(restricted.data$range.size.m2)
pairs(restricted.data[,c("log.range.size", "bioclim4", "residuals.PC1", "PC1.LIG", "NPP")], lower.panel=panel.smooth, upper.panel=panel.cor)

Figure S3: The highest correlation is between Long term temperature variation and seasonal variation. However this correlation is not excessively high (0.53) so we can be satisfied in including these predictors within our model.

Sexual Dimorphism

The sexual dimorphism dataset looks to have an overdispersed distribution (see phylogenic plot in manuscript). Unfortunately, transformations do not greatly improve the distribution. Our model fit is expected to be reduced by this but further transformations are not obvious and would risk problems of interperetation.


Analysis

Simple speciation measures: Diversification rate (DR) and Node Density (ND)

Diversification rates can be simply estimated from thew ‘tippiness of the phylogeny’. Here a log equal splits (also called DR) value can be achieved that is weighted to more recent diversification. It was used by Jetz et al. (2012) and recently (in comparison with BAMM) in a study ivestigation bird diversity alongside elevation (Quintero and Jetz 2018). We can look at speciation rate using ES-Sim adapted by Harvey Michael et al. (2017) and found here https://github.com/mgharvey/ES-sim.

The functions that generate these tip rates estimates are: calculate.log.es, calculate.nd and calculate.tb. Here they are run on 100 trees plus a full MCC tree of passerine birds based on 2500 samples of the posterior (Jetz et al. 2012). This was done using the phanghorn (Schliep 2010).

#Obtain DR (log(es)) estimates, by calling the calculate.dr function
passerine.trees <- read.nexus('bamm_extinction/data/trees/Passerine_Trees_Full.nex')

#DR, ND and TB can be calculated on the 100 trees by: 
  
# es.list <- sapply(passerine.trees, calculate.log.es)
# nd.list <- sapply(passerine.trees, calculate.nd)
# tb.list <- sapply(passerine.trees, calculate.tb)

#To avoid re-running the time-consuming rates, we can save them and load them: 
es.list <- readRDS("data/es.list.rds")
nd.list <- readRDS("data/nd.list.rds")
tb.list <- readRDS("data/tb.list.rds")

#Also load in our MCC tree
MCC.passerine <- read.tree('data/MCC.passerine.tre') #Tree based on 2500 samples of posterior

To gain an estimate of variation between trees we can take the values and obtain a mean, variance and coeeficient of variation (CV). We obtain a CV value to account for unequal variation at higher estimates. The CV is expressed as the ratio of the standard deviation (\(\sigma\)) to the mean (\(\mu\)) (converting to a percentage not needed):

\(C_V = \frac{\sigma}{\mu} \times 100\)

Furthermore, we can account for overdispersion by obtaining a transformed version of CV. In such a case we can calculate CV for a log-normal distribution:

\({\widehat {c_{\rm {v}}}}_{\rm {ln}}={\sqrt {\mathrm {e} ^{s_{\rm {ln}}^{2}}-1}}\)

Given that our ES value is log-transformed, the calculation of \({\widehat {c_{\rm {v}}}}_{\rm {ln}}\) is perhaps more ideal. As for ND we can see from the a qqnorm plot that value of ND is normally distributed, therefore the regular CV is more appropriate for this distribution.

#For ND
#transpose the list so that the elements are the species
t.nd.list <- nd.list %>% purrr::transpose()
#not too sure why this works but turn the nested list into a dataframe
nd.values <- as.data.frame(t(sapply(t.nd.list, function(x) sapply(x, function(x) (x)))))
nd.values <- nd.values %>% tibble::rownames_to_column("binomial")
nd.summary <- melt(nd.values, id.vars = "binomial") %>% group_by(binomial) %>% summarise(
  mean.nd = mean(value), #use normal mean, not a rate
  var.nd = var(value), #Don't really need this if we use CV
  CV.nd = (sd(value)/mean(value)) #Given normal distribution of ND, standard CV formula is appropriate
)

#For ES
#transpose the list so that the elements are the species
t.es.list <- es.list %>% purrr::transpose()
#not too sure why this works but turn the nested list into a dataframe
es.values <- as.data.frame(t(sapply(t.es.list, function(x) sapply(x, function(x) (x)))))
es.values <- es.values %>% tibble::rownames_to_column("binomial")
es.summary <- melt(es.values, id.vars = "binomial") %>% group_by(binomial) %>% summarise(
  mean.loges = log(1/mean(1/exp(value))), #use harmonic mean for rates and log-transform post calculation of HM
  var.loges = var(value), #Don't really need this if we use CV
  CV.loges = sqrt(exp(var(value)) - 1) #CV for log-normal distribution of ES
)

#For TB
#transpose the list so that the elements are the species
t.tb.list <- tb.list %>% purrr::transpose()
#not too sure why this works but turn the nested list into a dataframe
tb.values <- as.data.frame(t(sapply(t.tb.list, function(x) sapply(x, function(x) (x)))))
tb.values <- tb.values %>% tibble::rownames_to_column("binomial")
tb.summary <- melt(tb.values, id.vars = "binomial") %>% group_by(binomial) %>% summarise(
  mean.tb = mean(value), #use normal mean; not a rate
  var.tb = var(value), #Don't really need this if we use CV
  CV.tb = (sd(value)/mean(value)) #Given normal distribution of tb, standard CV formula is appropriate
)

MCC.nd.df <- as.data.frame(calculate.nd(MCC.passerine)) %>% tibble::rownames_to_column("binomial")
MCC.dr.df <- as.data.frame(calculate.log.es(MCC.passerine)) %>% tibble::rownames_to_column("binomial")

#Join the dataframes
dr.summary <- full_join(nd.summary, es.summary, by = "binomial")
dr.summary <- full_join(dr.summary, tb.summary, by = "binomial")
dr.summary <- full_join(dr.summary, MCC.nd.df, by = "binomial")
dr.summary <- full_join(dr.summary, MCC.dr.df, by = "binomial")

dr.summary <- dr.summary %>% rename(TipLabel = binomial, MCC.ND = nd, MCC.DR = es)

saveRDS(dr.summary, 'data/dr.summary.rds')
dr.summary <- readRDS('data/dr.summary.rds')
#Plot a summary of logES and TB with the DRs with their weightings
dr.plot <- dr.summary %>% ggplot(aes(x = mean.loges, y = mean.nd, 
                          fill = (1/(CV.nd))/sum((1/(CV.nd)))*100, 
                      size = (1/(CV.loges))/sum((1/(CV.loges)))*100))+
  geom_point(shape = 21, alpha = 0.5)+
  theme_minimal()+
  geom_errorbarh(aes(xmin = mean.loges - 0.674*sqrt(var.loges), xmax = mean.loges + 0.674*sqrt(var.loges)), 
                 size = 0.015)+
  geom_errorbar(aes(ymin = mean.nd - 0.674*sqrt(var.nd), ymax = mean.nd + 0.674*sqrt(var.nd)), 
                size = 0.015)+
  theme(legend.position="bottom")+
  labs(size = 'Weight (%) [using DR]', y='Node Density [ND]', x= 'Log Diversification Rate [DR]', fill = 'Weight (%) [using ND] \n')+
  scale_fill_distiller(guide = "colorbar", palette = "RdPu", direction = 1, breaks=c(0,0.015,.03),
                           limits=c(0,0.03))
marginal.dr <-ggMarginal(
  p = dr.plot,
  type = 'density',
  margins = 'both',
  size = 5,
  colour = '#0000009C',
  fill = '#D12E6769'
)
grid.newpage()
grid.draw(marginal.dr)

Figure S4: \(\lambda_{DR}\) and \(\lambda_{ND}\) are correlated in their estimation of species rate. Furthermore, we see similar variation across each species in the estimates across 100 trees as the inverse of the coefficeint of variation (weight) is approximately equivalent across all values of ND and DR (colours and size). Notably, these estimates speciation rate are relatively uncertain across all species, with 50 % CIs in both axes plotted.

Combine the estimates of DR and ND with the restricted dataframe containing all the other variables. We can inspect the trend in SDi in relation to the speciation rates.

restricted.data <- left_join(restricted.data, dr.summary, by = 'TipLabel')

es.plot <- restricted.data %>% ggplot(aes(y = mean.loges, x = SDi, fill = mean.loges))+
geom_point(shape = 21, size = 1.5)+
geom_smooth(show.legend = FALSE, color = "grey20", method = "lm")+
scale_fill_distiller(palette = "YlOrBr", direction = 1, guide = FALSE)+
ylab("mean.DR")+
theme_minimal()

nd.plot <- restricted.data %>% ggplot(aes(y = mean.nd, x = SDi, fill = mean.nd))+
geom_point(shape = 21, size = 1.5)+
geom_smooth(show.legend = FALSE, color = "grey20", method = "lm")+
scale_fill_distiller(palette = "Greens", direction = 1, guide = FALSE)+
theme_minimal()

tb.plot <- restricted.data %>% ggplot(aes(y = mean.tb, x = SDi, fill = mean.tb))+
geom_point(shape = 21, size = 1.5)+
geom_smooth(show.legend = FALSE, color = "grey20", method = "lm")+
scale_fill_distiller(palette = "PuBu", direction = 1, guide = FALSE)+
theme_minimal()

grid.arrange(es.plot, nd.plot, tb.plot, nrow = 3)

Figure S5: Scatter plots showing the raw relationship between sexual dichromatism (SDi) and speciation rates. Across all three measures the pattern and spread is similar, with no obvious relationship. Note that terminal branch lenght (mean.tb) was not used in the analysis.

BAMM measures of speciation and extinction

Set up BAMM parameters

For the use of BAMM we used the following code to generate the parameters across the 100 trees. Each paramaeter value is specified in this code chunk to ensure reproducibility when re-running on any of the trees. These same parameters were also used for the MCC tree (with the seed set at 2500).

name.passerine.tree <- names(passerine.trees)

priors <- sapply(name.passerine.tree, function(x) {
  setBAMMpriors(passerine.trees[[x]], outfile = NULL)
})

sapply(name.passerine.tree, function(x) {
  write.tree(passerine.trees[[x]], paste("data/bamm_files/", x, ".tre", sep=""))
})

# Here is a block of parameters for the control file. We can make a control file for each tree:
params <- list()
for (x in name.passerine.tree) {

# GENERAL SETUP AND DATA INPUT

params[[x]] <- list(modeltype = 'speciationextinction',
# Specify "speciationextinction" or "trait" analysis
                                  
treefile = paste(x, ".tre", sep=""),
# File name of the phylogenetic tree to be analyzed

runInfoFilename = 'run_info.txt',
# File name to output general information about this run

sampleFromPriorOnly = 0,
# Whether to perform analysis sampling from prior only (no likelihoods computed)

runMCMC = 1,
# Whether to perform the MCMC simulation. If runMCMC = 0, the program will only
# check whether the data file can be read and the initial likelihood computed

loadEventData = 0,                       
# Whether to load a previous event data file

eventDataInfile = 'event_data_in.txt',
# File name of the event data file to load, used only if loadEventData = 1

initializeModel = 1,
# Whether to initialize (but not run) the MCMC. If initializeModel = 0, the
# program will only ensure that the data files (e.g., treefile) can be read

useGlobalSamplingProbability = 1,
# Whether to use a "global" sampling probability. If False (0), expects a file
# name for species-specific sampling probabilities (see sampleProbsFilename)
                                        
globalSamplingFraction = 1,
# The sampling probability. If useGlobalSamplingProbability = 0, this is ignored
# and BAMM looks for a file name with species-specific sampling fractions

sampleProbsFilename = 'sample_probs.txt',
# File name containing species-specific sampling fractions

seed = as.numeric(gsub("tree_", "", x, perl = TRUE)),
# Seed for the random number generator. Set for reproducibility to the number of the treefile

overwrite = 1,
# If True (1), the program will overwrite any output files in the current
# directory (if present)


# PRIORS

expectedNumberOfShifts = 100,
# prior on the number of shifts in diversification
# Suggested values: 
#     expectedNumberOfShifts = 1.0 for small trees (< 500 tips)
#  expectedNumberOfShifts = 10 or even 50 for large trees (> 5000 tips) 
 
lambdaInitPrior = as.numeric(priors['lambdaInitPrior', x]),
# Prior (rate parameter of exponential) on the initial lambda value for rate
# regimes

lambdaShiftPrior = 0.05,
# Prior (std dev of normal) on lambda shift parameter for rate regimes
# You cannot adjust the mean of this distribution (fixed at zero, which is
# equal to a constant rate diversification process)

muInitPrior = as.numeric(priors['muInitPrior', x]),
# Prior (rate parameter of exponential) on extinction rates  

lambdaIsTimeVariablePrior = 1,
# Prior (probability) of the time mode being time-variable (vs. time-constant)
            

# MCMC SIMULATION SETTINGS & OUTPUT OPTIONS

numberOfGenerations = '100000000',
# Number of generations to perform MCMC simulation

mcmcOutfile = 'mcmc_out.txt',
# File name for the MCMC output, which only includes summary information about
# MCMC simulation (e.g., log-likelihoods, log-prior, number of processes)

mcmcWriteFreq = 1000,
# Frequency in which to write the MCMC output to a file

eventDataOutfile = 'event_data.txt',
# The raw event data (these are the main results). ALL of the results are
# contained in this file, and all branch-specific speciation rates, shift
# positions, marginal distributions etc can be reconstructed from this output.
# See R package BAMMtools for working with this output

eventDataWriteFreq = 1000,
# Frequency in which to write the event data to a file

printFreq = 10000,
# Frequency in which to print MCMC status to the screen

acceptanceResetFreq = 1000,
# Frequency in which to reset the acceptance rate calculation
# The acceptance rate is output to both the MCMC data file and the screen

outName = x,
# Optional name that will be prefixed on all output files (separated with "_")
# If commented out, no prefix will be used


# OPERATORS: MCMC SCALING OPERATORS

updateLambdaInitScale = 2,
# Scale parameter for updating the initial speciation rate for each process

updateLambdaShiftScale = 0.1,
# Scale parameter for the exponential change parameter for speciation

updateMuInitScale = 2,
# Scale parameter for updating initial extinction rate for each process

updateEventLocationScale = 0.1,
# Scale parameter for updating LOCAL moves of events on the tree
# This defines the width of the sliding window proposal
 
updateEventRateScale = 4,
# Scale parameter (proportional shrinking/expanding) for updating
# the rate parameter of the Poisson process

# OPERATORS: MCMC MOVE FREQUENCIES

updateRateEventNumber = 1,
# Relative frequency of MCMC moves that change the number of events

updateRateEventPosition = 0.25,
# Relative frequency of MCMC moves that change the location of an event on the
# tree

updateRateEventRate = 1,
# Relative frequency of MCMC moves that change the rate at which events occur 

updateRateLambda0 = 1,
# Relative frequency of MCMC moves that change the initial speciation rate
# associated with an event

updateRateLambdaShift = 1,
# Relative frequency of MCMC moves that change the exponential shift parameter
# of the speciation rate associated with an event

updateRateMu0 = 1,
# Relative frequency of MCMC moves that change the extinction rate for a given
# event

updateRateLambdaTimeMode = 0,
# Relative frequency of MCMC moves that flip the time mode
# (time-constant <=> time-variable)

localGlobalMoveRatio = 10,
# Ratio of local to global moves of events 


# INITIAL PARAMETER VALUES

lambdaInit0 = 0.032,
# Initial speciation rate (at the root of the tree)

lambdaShift0 = 0,
# Initial shift parameter for the root process

muInit0 = 0.005,
# Initial value of extinction (at the root)

initialNumberEvents = 0,
# Initial number of non-root processes


# METROPOLIS COUPLED MCMC

numberOfChains = 1,
# Number of Markov chains to run

deltaT = 0.01,
# Temperature increment parameter. This value should be > 0
# The temperature for the i-th chain is computed as 1 / [1 + deltaT * (i - 1)]

swapPeriod = 1000,
# Number of generations in which to propose a chain swap

chainSwapFileName = 'chain_swap.txt',
# File name in which to output data about each chain swap proposal.
# The format of each line is [generation],[rank_1],[rank_2],[swap_accepted]
# where [generation] is the generation in which the swap proposal was made,
# [rank_1] and [rank_2] are the chains that were chosen, and [swap_accepted] is
# whether the swap was made. The cold chain has a rank of 1.


# NUMERICAL AND OTHER PARAMETERS

minCladeSizeForShift = 3,
# Allows you to constrain location of possible rate-change events to occur
# only on branches with at least this many descendant tips. A value of 1
# allows shifts to occur on all branches. 

segLength = 0.025,
# Controls the "grain" of the likelihood calculations. Approximates the
# continuous-time change in diversification rates by breaking each branch into
# a constant-rate diversification segments, with each segment given a length
# determined by segLength. segLength is in units of the root-to-tip distance of
# the tree. So, if the segLength parameter is 0.01, and the crown age of your
# tree is 50, the "step size" of the constant rate approximation will be 0.5.
# If the value is greater than the branch length (e.g., you have a branch of
# length < 0.5 in the preceding example) BAMM will not break the branch into
# segments but use the mean rate across the entire branch.

outName = x)
  }

bammcontrolfile <- list()
for (x in name.passerine.tree) {
  bammcontrolfile[x] <- paste("data/bamm_files/control_", x, ".txt", sep="")
}

# Now writing control parameters to file

for (x in name.passerine.tree) {generateControlFile(file = bammcontrolfile[[x]], type = "diversification", params = params[[x]])}

Run analysis

BAMM can be run through the terminal through the following syntax: bamm -c control_tree_xxxx.txt. To generate these commands we can use a loop function, from which we get:

bamm.commands <- list()
for (x in name.passerine.tree) {
  bamm.commands[x] <- paste("bamm -c control_", x, ".txt", sep="")
}

The analysis was run over multiple CPU’s, each generating a respective MCMC and EventData output. Due to the size of the event data file (~ 50 Gb in total) they are not included as supplementary material here. However we have simplified the event data objects into tip rate estimates of the mean and variance across 100 trees + MCC.

Read in the event data and extract tip data

Now we have a series of trees and event data we can read in. Firstly and most importantly let’s check for convergence in the MCC tree through the following:

Table S3: Effective sample size (ESS) for the two key BAMM parameters (number of evolutionary shifts and log-Likelihood) for the run on the MCC indicate that BAMM converges (ESS > 200).

#Read in the tree and MCMC to check for convergence 
MCC.BAMM.tree  <- read.tree("data/BAMM_MCC/MCC.passerine.tre") #Same as other MCC tree already loaded
MCC.BAMM.mcmc <- read.csv( "data/BAMM_MCC/tree_MCC_mcmc_out.txt" , stringsAsFactors=F)

#Plot of convergence can be generated by:
#plot(MCC.BAMM.mcmc$logLik ~ MCC.BAMM.mcmc$generation)

#Looks like it has converged so let's discard burn in: 
burnstart <- floor(0.1 * nrow(MCC.BAMM.mcmc))
postburn <- MCC.BAMM.mcmc[burnstart:nrow(MCC.BAMM.mcmc), ]

#We can also check effective population sizes of the log-likelihhod and number of shift events in each sample
#We want at least 200 (although that's on the low side)

cbind(effectiveSize(postburn$N_shifts), effectiveSize(postburn$logLik)) %>% `colnames<-`(c("N_Shifts", "logLik")) %>% `rownames<-`("Effective Sample Size") %>% pander()
  N_Shifts logLik
Effective Sample Size 2456 1302

Those values are well above 200, so we are good to use the MCC data in analysis! We can also check the convergence of BAMM across 100 runs of BAMM. The Raw MCMCs are not included in this file or repository but we can read in a dataframe that has extracted the ESS for each of the runs.

Table S4: For the 100 trees that BAMM was run on effective sample size (ESS) for the two key BAMM parameters (number of evolutionary shifts and log-Likelihood) also indicates that BAMM converges with the minimum for each parameter across the 100 trees being over 200.

ESS <- readRDS('data/ESS.rds')
summary(ESS) %>% pander()
N_Shifts logLik
Min. : 723.4 Min. : 278.7
1st Qu.:1396.1 1st Qu.: 511.3
Median :1618.6 Median : 638.6
Mean :1709.5 Mean : 669.9
3rd Qu.:2070.8 3rd Qu.: 804.8
Max. :3425.6 Max. :1244.1

Given that it appears BAMM converges we can we can make a data frame with the mean and variance for extinction and speciation tip rates from the large event data set for the MCC with the following code and then plot the variation we see in tip-rates.

# Read in Event Data
MCC.BAMM.ED <- getEventData(MCC.BAMM.tree,  "data/BAMM_MCC/tree_MCC_event_data.txt", burnin=0.1, nsamples=1000)
## Reading event datafile:  data/BAMM_MCC/tree_MCC_event_data.txt 
##      ...........
## Read a total of 100000 samples from posterior
## 
## Discarded as burnin: GENERATIONS <  9999000
## Analyzing  1000  samples from posterior
## 
## Setting recursive sequence on tree...
## 
## Done with recursive sequence
saveRDS(MCC.BAMM.ED, 'data/MCC.BAMM.ED.rds')
#From the Event Data we can extract 

library(purrr)
#BAMM.EventData <- readRDS('data/BAMM.EventData.rds')

#Big lapply over each tree in BAMM event data
BAMM.extraction.function <- function(x) {
######Get mean and var for lambda
#Transpose list so each element in the list is a species
transposed.lambda <- lapply(purrr::transpose(x$tipLambda), unlist)

#Now turn it into a df with mean and variance
lambda <- sapply(transposed.lambda, function(x) {
  mean.lambda = mean(log(x))
  var.lambda = var(log(x))
  return(c(mean.lambda, var.lambda))
}) %>% t() %>% as.data.frame() %>% `colnames<-`(c("mean.lambda", "var.lambda"))

lambda$TipLabel <- x[["tip.label"]]

#####NOW FOR Extinction

#Transpose list so each element in the list is a species
transposed.mu <- lapply(purrr::transpose(x$tipMu), unlist)

#Now turn it into a df with mean and variance
mu <- sapply(transposed.mu, function(x) {
  mean.mu = mean(log(x))
  var.mu = var(log(x))
  return(c(mean.mu, var.mu))
}) %>% t() %>% as.data.frame() %>% `colnames<-`(c("mean.mu", "var.mu"))

mu$TipLabel <- x[["tip.label"]]

left_join(lambda, mu, by = "TipLabel")
}

MCC.BAMM.df <- BAMM.extraction.function(MCC.BAMM.ED)

#Save df for later use
saveRDS(MCC.BAMM.df, 'data/MCC.BAMM.df.rds')
MCC.BAMM.df <- readRDS('data/MCC.BAMM.df.rds')
#Plot a summary of logES and TB with the DRs with their weightings
BAMM.MCC.plot <- MCC.BAMM.df %>% ggplot(aes(x = mean.lambda, y = mean.mu, 
                          fill = 1/var.lambda, 
                      size = 1/var.mu))+
  geom_point(shape = 21, alpha = 0.5)+
  theme_minimal()+
  geom_errorbarh(aes(xmin = mean.lambda - 0.674*sqrt(var.lambda), xmax = mean.lambda + 0.674*sqrt(var.lambda)),
                 size = 0.0025)+
  geom_errorbar(aes(ymin = mean.mu - 0.674*sqrt(var.mu), ymax = mean.mu + 0.674*sqrt(var.mu)),
                size = 0.0025)+
  # scale_y_continuous(trans = "log")+
  # scale_x_continuous(trans = "log")+
  # xlim(-8, 3)+
  # ylim(-10,2)+
  # scale_x_continuous(trans = "log")+
  # scale_y_continuous(trans = "log")+
  theme(legend.position="bottom")+
  labs(size = 'Inverse log(var) / Weight [using Lambda]', y='Log Extinction [Mu]', x= 'Log Speciation [Lambda]', fill = 'Inverse log(var) / Weight [using Mu] \n')+
  scale_fill_distiller(guide = "colorbar", palette = "Reds", direction = 1)

BAMM.variance <- ggExtra::ggMarginal(
  p = BAMM.MCC.plot,
  type = 'density',
  margins = 'both',
  size = 5,
  colour = 'black',
  fill = '#BA3B1C91'
)
grid.newpage()
grid.draw(BAMM.variance)

Figure S5: The tip-rate estimate for BAMM are highly variable within each run of BAMM. Across most species there is high variability in the posterior distribution of tip-rate estimates. Here we show mean values, weights based on the variance and 50 % CIs.

Analysis of BAMM results

Given the high variability in tip-rate estimates from the above plot, below we peform some diagnostics on BAMM to demonstrate that the variability is unlikely an error in sampling or parameters, rather an inherent aspect of BAMM. The following is not inherently necessary to understand the conclusions drawn from the paper, however they do raise a set of methodological questions about BAMM that warrant further investigation.

Credible number of shifts

To plot the credible shift set, we need the prior distribution on the number of rate shifts (this is generated internally by BAMMtools). We can then estimate the credible set of rate shifts using the BAMMtools function credibleShiftSet:

css <- credibleShiftSet(MCC.BAMM.ED, expectedNumberOfShifts=100, threshold=5, set.limit = 0.95)

#Now we obtain the number of distinct shifts: (Out of 1000 samples this is super high, essentially each one is distinct)

css$number.distinct
## [1] 950

This value indicates that from a posterior sample of 1000 we have 950 unique rate shifts, which is very high. Arguably, across such a large phylogenetic tree, even if the number of shifts reaches convergence the location of those > 50 shifts will be quite variable, with the variation likely impacting tip-rate estimates. Even if we increase the sampling of the posterior to 2,500 instead of 1,000 we still see that the CSS is high (close to 2500). We can have a look what our top few shifts are:

summary(css)
## 
##  95 % credible set of rate shift configurations sampled with BAMM
## 
## Distinct shift configurations in credible set:  950
## 
## Frequency of 9 shift configurations with highest posterior probability:
## 
## 
##    rank     probability cumulative  Core_shifts
##          1      0.001      0.001         54
##          2      0.001      0.002         44
##          3      0.001      0.003         46
##          4      0.001      0.004         47
##          5      0.001      0.005         49
##          6      0.001      0.006         45
##          7      0.001      0.007         47
##          8      0.001      0.008         43
##          9      0.001      0.009         45
## 
## ...omitted 941 additional distinct shift configurations
## from the credible set. You can access the full set from your 
## credibleshiftset object

Notably, each of the top shifts has low probability, indicating that out of the entire posterior sample we cannot differentiate which shift configuration is more likely. Based on the follwing eastimates of the BayesFactor, we are convinced that the number of shifts is non-zero, however it is still quite a wide distribution, sourced from potentially weakly non-informative priors. But, this is at odds with the ESS and the following plot:

round(computeBayesFactors(MCC.BAMM.mcmc, expectedNumberOfShifts=100, burnin=0.1)[,1], digits = 2)
##      42      43      44      45      46      47      48      49      50 
##    1.00   17.17   24.48   49.45  120.71  283.77  454.33  784.80 1291.85 
##      51      52      53      54      55      56      57      58      59 
## 1813.33 2561.62 3608.07 4574.91 5492.44 6322.11 7141.12 7720.26 7886.28 
##      60      61      62      63      64      65      66      67      68 
## 7846.73 7661.83 7116.15 6466.36 5749.34 4811.16 4136.80 3338.17 2674.70 
##      69      70      71      72      73      74      75      76      77 
## 2137.61 1593.48 1218.40  868.01  643.91  452.36  311.07  199.17  135.99 
##      78      79      80      81      82      83      84      85      87 
##   85.85   65.03   40.87   14.74   13.40    9.02   10.63    1.53    1.56
plotPrior(MCC.BAMM.mcmc, expectedNumberOfShifts=100)

Figure S6: The apparant convergence of the number of shifts in the posterior is at odds with the variability seen in the CSS and although there is greater certainty in the number of shifts being within the range above, the position of those shifts remain highly variable.

We can compare our run of BAMM against Harvey et al. (2017) who used BAMM on a genetic-only MCC tree with different parameters

load('data/Hackett_split_eventsample.rda')
css3 <- credibleShiftSet(ed, expectedNumberOfShifts=100, threshold=5, set.limit = 0.95)
css3$number.distinct
## [1] 1188

Based on Their posterior sample of 1250, 1188 are unique. Like our run on the MCC, we find high uncertainity in the position of shift-configurations.

Let’s run a model based on genetic data from the Harvey MCC

Harvey.BAMM.df <- BAMM.extraction.function(ed)

Harvey.BAMM.df %>% ggplot(aes(x = mean.lambda, y = mean.mu, 
                          fill = 1/(var.lambda), 
                      size = 1/(var.mu)))+
  geom_point(shape = 21, alpha = 0.5)+
  theme_minimal()+
  geom_errorbarh(aes(xmin = mean.lambda - 0.674*sqrt(var.lambda), xmax = mean.lambda + 0.674*sqrt(var.lambda)), 
                 size = 0.0025)+
  geom_errorbar(aes(ymin = mean.mu - 0.674*sqrt(var.mu), ymax = mean.mu + 0.674*sqrt(var.mu)), 
                size = 0.0025)+
  # scale_x_continuous(trans = "log")+
  # scale_y_continuous(trans = "log")+
  theme(legend.position="bottom")+
  labs(size = 'Weight [using Lambda]', y='Log Extinction [Mu]', x= 'Log Speciation [Lambda]', fill = 'Weight [using Mu] \n')+
  scale_fill_distiller(guide = "colorbar", palette = "Reds", direction = 1)

# BAMM.variance <- ggExtra::ggMarginal(
#   p = BAMM.MCC.plot,
#   type = 'density',
#   margins = 'both',
#   size = 5,
#   colour = 'black',
#   fill = '#BA3B1C91'
# )

Figure S7: Based on the event data from the BAMM run by Harvey et al. (2017) we find that in both our case and theirs the tip-rate estimates for many species are extrememly variable. here we present mean estimates with 50 % CIs.

Time-constant

We also ran a time-constant MCC of BAMM

#Read in the tree and MCMC to check for convergence 
MCC.BAMM.tree  <- read.tree("data/BAMM_MCC/MCC.passerine.tre") #Same as other MCC tree already loaded
MCC.BAMM.Constant.mcmc <- read.csv( "data/BAMM_MCC/TimeConstant/tree_MCC_Constant_mcmc_out.txt" , stringsAsFactors=F)

#Plot of convergence can be generated by:
#plot(MCC.BAMM.mcmc$logLik ~ MCC.BAMM.mcmc$generation)

#Looks like it has converged so let's discard burn in: 
burnstart <- floor(0.1 * nrow(MCC.BAMM.Constant.mcmc))
postburn <- MCC.BAMM.Constant.mcmc[burnstart:nrow(MCC.BAMM.Constant.mcmc), ]

#We can also check effective population sizes of the log-likelihhod and number of shift events in each sample
#We want at least 200 (although that's on the low side)

cbind(effectiveSize(postburn$N_shifts), effectiveSize(postburn$logLik)) %>% `colnames<-`(c("N_Shifts", "logLik")) %>% `rownames<-`("Effective Sample Size") %>% pander()


# Read in Event Data
MCC.BAMM.ED.Constant <- getEventData(MCC.BAMM.tree,  "data/BAMM_MCC/TimeConstant/tree_MCC_Constant_event_data.txt", burnin=0.1, nsamples=1000)
saveRDS(MCC.BAMM.ED.Constant, 'data/MCC.BAMM.ED.Constant.rds')
#From the Event Data we can extract 

css2 <- credibleShiftSet(MCC.BAMM.ED.Constant, expectedNumberOfShifts=100, threshold=5, set.limit = 0.95)

#Now we obtain the number of distinct shifts: (Out of 1000 samples this is super high, essentially each one is distinct)

css2$number.distinct

plotPrior(MCC.BAMM.Constant.mcmc, expectedNumberOfShifts=100)

MCC.BAMM.df.Constant <- BAMM.extraction.function(MCC.BAMM.ED.Constant)

MCC.BAMM.df.Constant %>% ggplot(aes(x = mean.lambda, y = MCC.BAMM.df$mean.lambda))+
  geom_point()
cor(MCC.BAMM.df.Constant$mean.lambda, MCC.BAMM.df$mean.lambda)

PGLS Models

The method behind running the PGLS models is as follows:

1: Estimate the phylogenetic signal (\(\lambda\)) by running a model without interactions and all six predictor variables. The value of \(\lambda\) obtained here will be fixed in all successive models. The value is fixed for subsequent models as independently estimating it in each case becomes computationally intensive. Based on preliminary analysis even slight differences in lambda lead to large reductions in AIC.

2: Create a global model of six predictor variables plus the fiver interactions between sexual dichromatism and the other variables.

3: Dredge the global model but fix the six independent predictors, hence conducting model selection on the interaction terms.

4: Take the top model in the MCC model and run it on the 100 phylogenetic trees.

5: Repeat this step for DR, ND, BAMM-speciation and BAMM-extinction.

PGLS Models on DR and ND

Using corPagel we can estimate the phylogenetic signal for a model with all predictors (interactions do not appear to affect the estimate of \(\lambda\)):

#Prune tree
pruned.MCC.tree <- drop.tip(MCC.passerine,MCC.passerine$tip.label[-match(restricted.data$TipLabel, MCC.passerine$tip.label)])

#Set rownames to match tree
rownames(restricted.data) <- restricted.data$TipLabel

#Run a corPagel model to estimate lambda for DR
MCC.DR.corPagel <- gls(MCC.DR ~ SDi
                         + log(range.size.m2)
                         + bioclim4 #Seasonal variation
                         + residuals.PC1 #Spatial variation
                         + PC1.LIG #Long-term climate variation
                         + NPP,
                correlation = corPagel(1, phy = pruned.MCC.tree, fixed = FALSE), 
                data = restricted.data, 
                method = "REML")
saveRDS(MCC.DR.corPagel, 'data/MCC.DR.corPagel.rds')

#Run a corPagel model to estimate lambda for ND
MCC.ND.corPagel <- gls(MCC.ND ~ SDi
                         + log(range.size.m2)
                         + bioclim4 #Seasonal variation
                         + residuals.PC1 #Spatial variation
                         + PC1.LIG #Long-term climate variation
                         + NPP,
                correlation = corPagel(1, phy = pruned.MCC.tree, fixed = FALSE), 
                data = restricted.data, 
                method = "REML")
saveRDS(MCC.ND.corPagel, 'data/MCC.ND.corPagel.rds')

Inspect the \(\lambda\) value, which can then be fixed for successive models.

MCC.DR.corPagel <- readRDS('data/MCC.DR.corPagel.rds')
MCC.ND.corPagel <- readRDS('data/MCC.ND.corPagel.rds')
MCC.DR.corPagel[["modelStruct"]][["corStruct"]] %>% `names<-`("DR lambda") %>% pander()
  • DR lambda: 0.985
MCC.ND.corPagel[["modelStruct"]][["corStruct"]] %>% `names<-`("ND lambda") %>% pander()
  • ND lambda: 0.9996

The lambda is high and not too dissimilar than if we assume Brownian Motion \(\lambda = 1\). However, given the large sample size this difference may have an effect on the results so, to play safe we included it as a fixed value for \(\lambda\) in all successive models.

#Run model for DR
MCC.DR <- gls(MCC.DR ~ SDi
                         + log(range.size.m2)
                         + bioclim4 #Seasonal variation
                         + residuals.PC1 #Spatial variation
                         + PC1.LIG #Long-term climate variation
                         + NPP
                         + SDi*log(range.size.m2)
                         + SDi*bioclim4
                         + SDi*residuals.PC1
                         + SDi*PC1.LIG
                         + SDi*NPP,
                correlation = corPagel(0.985, phy = pruned.MCC.tree, fixed = TRUE), 
                data = restricted.data, 
                method = "REML")

#Run model for ND
MCC.ND <- gls(MCC.ND ~ SDi
                      + log(range.size.m2)
                      + bioclim4 #Seasonal variation
                      + residuals.PC1 #Spatial variation
                      + PC1.LIG #Long-term climate variation
                      + NPP
                      + SDi*log(range.size.m2)
                      + SDi*bioclim4
                      + SDi*residuals.PC1
                      + SDi*PC1.LIG
                      + SDi*NPP,
                correlation = corPagel(0.9996, phy = pruned.MCC.tree, fixed = TRUE), 
                data = restricted.data, 
                method = "REML")

#Set up cluster
cores<-8
clusterType <- if(length(find.package("snow", quiet = TRUE))) "SOCK" else "PSOCK"
clust <- try(makeCluster(getOption("cl.cores", cores), type = clusterType))
myclust<-clust

#Export data and packages to cluster
clusterExport(myclust, c("restricted.data"), envir=environment())
clusterExport(myclust, c("pruned.MCC.tree"), envir=environment())
clusterEvalQ(myclust, library(nlme))
clusterEvalQ(myclust, library(ape))
clusterEvalQ(myclust, library(MuMIn))

#Dredged models:

dredged.ND.model <- pdredge(MCC.ND, fixed = c("SDi", "log(range.size.m2)", "bioclim4", "residuals.PC1", "PC1.LIG", "NPP"), cluster=myclust, trace = TRUE)

saveRDS(dredged.ND.model, "data/dredged.ND.model.rds")

dredged.DR.model <- pdredge(MCC.DR, fixed = c("SDi", "log(range.size.m2)", "bioclim4", "residuals.PC1", "PC1.LIG", "NPP"), cluster=myclust, trace = TRUE)

saveRDS(dredged.DR.model, "data/dredged.DR.model.rds")

Table S5: The dredged models both show the top model is one with no interactions, with \(\delta AICc > 4\) in both cases. We can be resonably confident that interactions are unlikely to affect the pattern of speciation we see in passerine birds.

dredged.DR.model <- readRDS("data/dredged.DR.model.rds")
dredged.ND.model <- readRDS("data/dredged.ND.model.rds")
kable(dredged.DR.model, "html", caption = "DR Dredge Table") %>%
  kable_styling() %>%
  scroll_box(width = "100%", height = "500px")
DR Dredge Table
(Intercept) bioclim4 log(range.size.m2) NPP PC1.LIG residuals.PC1 SDi bioclim4:SDi log(range.size.m2):SDi NPP:SDi PC1.LIG:SDi residuals.PC1:SDi df logLik AICc delta weight
0 -2.879811 2.14e-05 -0.0065780 3e-07 0.0015864 -0.0024295 -0.0012787 NA NA NA NA NA 8 -3409.547 6835.118 0.00000 9.988652e-01
8 -2.879495 2.14e-05 -0.0065788 3e-07 -0.0001028 -0.0023765 -0.0012942 NA NA NA 0.0003487 NA 9 -3416.089 6850.208 15.09006 5.281325e-04
16 -2.879345 2.15e-05 -0.0066034 3e-07 0.0015509 -0.0015982 -0.0013082 NA NA NA NA -0.0001638 9 -3416.488 6851.007 15.88895 3.542138e-04
2 -2.884345 2.14e-05 -0.0064052 3e-07 0.0015991 -0.0024487 -0.0003035 NA -0.0000369 NA NA NA 9 -3416.842 6851.714 16.59633 2.486904e-04
1 -2.878432 1.57e-05 -0.0065688 3e-07 0.0015872 -0.0024180 -0.0015875 1.2e-06 NA NA NA NA 9 -3421.210 6860.452 25.33361 3.150533e-06
24 -2.878703 2.16e-05 -0.0066190 3e-07 -0.0004658 -0.0010549 -0.0013437 NA NA NA 0.0004121 -0.0002586 10 -3422.848 6865.734 30.61631 2.245227e-07
4 -2.884812 2.14e-05 -0.0065580 7e-07 0.0015560 -0.0024009 -0.0003587 NA NA -1e-07 NA NA 9 -3424.307 6866.644 31.52618 1.424567e-07
10 -2.890455 2.13e-05 -0.0061602 3e-07 -0.0001797 -0.0024196 0.0010665 NA -0.0000894 NA 0.0003710 NA 10 -3423.322 6866.681 31.56307 1.398533e-07
18 -2.885405 2.15e-05 -0.0063727 3e-07 0.0015667 -0.0015933 -0.0000022 NA -0.0000495 NA NA -0.0001699 10 -3423.771 6867.580 32.46202 8.922136e-08
9 -2.879477 2.13e-05 -0.0065786 3e-07 -0.0000996 -0.0023764 -0.0012983 0.0e+00 NA NA 0.0003480 NA 10 -3427.745 6875.529 40.41070 1.676618e-09
17 -2.876970 1.25e-05 -0.0065992 3e-07 0.0015377 -0.0012411 -0.0018098 1.9e-06 NA NA NA -0.0002306 10 -3428.008 6876.054 40.93602 1.289324e-09
3 -2.889434 1.34e-05 -0.0061257 3e-07 0.0016199 -0.0024624 0.0007743 1.7e-06 -0.0000939 NA NA NA 10 -3428.372 6876.782 41.66339 8.962234e-10
12 -2.884227 2.14e-05 -0.0065599 7e-07 -0.0000916 -0.0023507 -0.0004248 NA NA -1e-07 0.0003405 NA 10 -3430.871 6881.780 46.66167 7.362996e-11
26 -2.893478 2.15e-05 -0.0060556 3e-07 -0.0006022 -0.0009958 0.0018498 NA -0.0001211 NA 0.0004478 -0.0002816 11 -3430.028 6882.101 46.98312 6.269768e-11
20 -2.884117 2.15e-05 -0.0065817 7e-07 0.0015263 -0.0016648 -0.0004369 NA NA -1e-07 NA -0.0001454 10 -3431.269 6882.576 47.45753 4.945782e-11
6 -2.895008 2.14e-05 -0.0061863 7e-07 0.0015803 -0.0024391 0.0018191 NA -0.0000790 -1e-07 NA NA 10 -3431.545 6883.128 48.01036 3.751366e-11
25 -2.877827 1.81e-05 -0.0066164 3e-07 -0.0003444 -0.0009513 -0.0015348 7.0e-07 NA NA 0.0003863 -0.0002784 11 -3434.448 6890.942 55.82373 7.542887e-13
11 -2.891727 1.89e-05 -0.0060861 3e-07 -0.0000889 -0.0024252 0.0013324 5.0e-07 -0.0001044 NA 0.0003533 NA 11 -3434.894 6891.834 56.71575 4.828789e-13
5 -2.883424 1.60e-05 -0.0065495 7e-07 0.0015572 -0.0023903 -0.0006670 1.1e-06 NA -1e-07 NA NA 10 -3435.976 6891.989 56.87126 4.467550e-13
19 -2.894722 8.10e-06 -0.0058780 3e-07 0.0015805 -0.0010549 0.0020214 2.8e-06 -0.0001542 NA NA -0.0002814 11 -3435.056 6892.158 57.03960 4.106919e-13
28 -2.882980 2.16e-05 -0.0065992 6e-07 -0.0004288 -0.0011303 -0.0005650 NA NA -1e-07 0.0004000 -0.0002393 11 -3437.666 6897.377 62.25931 3.020435e-14
14 -2.901190 2.13e-05 -0.0059396 7e-07 -0.0002030 -0.0024099 0.0032037 NA -0.0001318 -1e-07 0.0003719 NA 11 -3438.023 6898.091 62.97260 2.114363e-14
22 -2.895506 2.14e-05 -0.0061663 7e-07 0.0015517 -0.0016637 0.0019999 NA -0.0000885 -1e-07 NA -0.0001541 11 -3438.494 6899.034 63.91561 1.319496e-14
27 -2.898432 1.32e-05 -0.0057816 3e-07 -0.0003798 -0.0007224 0.0029159 1.7e-06 -0.0001787 NA 0.0004037 -0.0003394 12 -3441.451 6906.955 71.83674 2.513946e-16
13 -2.884235 2.15e-05 -0.0065599 7e-07 -0.0000930 -0.0023508 -0.0004230 0.0e+00 NA -1e-07 0.0003408 NA 11 -3442.528 6907.101 71.98280 2.336900e-16
21 -2.881649 1.31e-05 -0.0065788 6e-07 0.0015150 -0.0013253 -0.0009513 1.8e-06 NA -1e-07 NA -0.0002092 11 -3442.806 6907.657 72.53841 1.770074e-16
7 -2.901278 1.25e-05 -0.0058625 8e-07 0.0016024 -0.0024539 0.0031392 1.8e-06 -0.0001449 -1e-07 NA NA 11 -3443.046 6908.138 73.01974 1.391463e-16
30 -2.903249 2.15e-05 -0.0058572 7e-07 -0.0005994 -0.0010691 0.0037845 NA -0.0001586 -1e-07 0.0004443 -0.0002653 12 -3444.762 6913.578 78.46031 9.163629e-18
29 -2.882138 1.85e-05 -0.0065971 6e-07 -0.0003202 -0.0010362 -0.0007469 7.0e-07 NA -1e-07 0.0003770 -0.0002574 12 -3449.269 6922.592 87.47390 1.011093e-19
15 -2.903204 1.79e-05 -0.0058308 7e-07 -0.0000759 -0.0024175 0.0036222 7.0e-07 -0.0001538 -1e-07 0.0003472 NA 12 -3449.583 6923.220 88.10200 7.385839e-20
23 -2.905869 7.40e-06 -0.0056354 7e-07 0.0015654 -0.0011030 0.0042394 2.9e-06 -0.0002007 -1e-07 NA -0.0002702 12 -3449.752 6923.557 88.43915 6.240085e-20
31 -2.909040 1.25e-05 -0.0055511 7e-07 -0.0003574 -0.0007748 0.0050239 1.9e-06 -0.0002227 -1e-07 0.0003961 -0.0003275 13 -3456.166 6938.395 103.27731 3.742176e-23
kable(dredged.ND.model, "html", caption = "ND Dredge Table") %>%
  kable_styling() %>%
  scroll_box(width = "100%", height = "500px")
ND Dredge Table
(Intercept) bioclim4 log(range.size.m2) NPP PC1.LIG residuals.PC1 SDi bioclim4:SDi log(range.size.m2):SDi NPP:SDi PC1.LIG:SDi residuals.PC1:SDi df logLik AICc delta weight
0 0.1404205 5e-07 -0.0001462 0 -0.0000102 -6.47e-05 -0.0000575 NA NA NA NA NA 8 14329.64 -28643.26 0.00000 9.999550e-01
8 0.1403564 4e-07 -0.0001422 0 -0.0000850 -5.71e-05 -0.0000618 NA NA NA 1.48e-05 NA 9 14320.05 -28622.06 21.19677 2.495513e-05
16 0.1404275 5e-07 -0.0001466 0 -0.0000108 -5.67e-05 -0.0000581 NA NA NA NA -1.6e-06 9 14319.25 -28620.46 22.79683 1.121276e-05
2 0.1400429 5e-07 -0.0001322 0 -0.0000107 -6.61e-05 0.0000205 NA -2.9e-06 NA NA NA 9 14318.99 -28619.95 23.30593 8.692847e-06
1 0.1403637 6e-07 -0.0001457 0 -0.0000101 -6.47e-05 -0.0000496 0e+00 NA NA NA NA 9 14314.58 -28611.13 32.12698 1.056075e-07
4 0.1406803 5e-07 -0.0001499 0 -0.0000090 -6.52e-05 -0.0000932 NA NA 0 NA NA 9 14311.61 -28605.19 38.06105 5.434109e-09
24 0.1403735 4e-07 -0.0001429 0 -0.0000910 -3.30e-05 -0.0000638 NA NA NA 1.57e-05 -4.8e-06 10 14309.74 -28599.43 43.82308 3.047333e-10
10 0.1398108 4e-07 -0.0001218 0 -0.0000886 -5.88e-05 0.0000503 NA -4.2e-06 NA 1.54e-05 NA 10 14309.46 -28598.88 44.37966 2.307068e-10
18 0.1400302 5e-07 -0.0001318 0 -0.0000116 -5.60e-05 0.0000242 NA -3.1e-06 NA NA -2.1e-06 10 14308.61 -28597.18 46.08010 9.858583e-11
9 0.1401845 8e-07 -0.0001400 0 -0.0000977 -5.59e-05 -0.0000401 -1e-07 NA NA 1.74e-05 NA 10 14305.25 -28590.46 52.79550 3.432282e-12
17 0.1403708 6e-07 -0.0001459 0 -0.0000104 -6.13e-05 -0.0000504 0e+00 NA NA NA -7.0e-07 10 14304.23 -28588.42 54.83371 1.238774e-12
3 0.1400821 5e-07 -0.0001346 0 -0.0000106 -6.59e-05 0.0000102 0e+00 -2.4e-06 NA NA NA 10 14304.00 -28587.96 55.29503 9.835981e-13
12 0.1406450 4e-07 -0.0001461 0 -0.0000867 -5.74e-05 -0.0001020 NA NA 0 1.54e-05 NA 10 14302.09 -28584.14 59.11661 1.455366e-13
20 0.1406914 5e-07 -0.0001503 0 -0.0000097 -5.54e-05 -0.0000943 NA NA 0 NA -2.0e-06 10 14301.22 -28582.41 60.84729 6.125736e-14
6 0.1404437 5e-07 -0.0001414 0 -0.0000094 -6.61e-05 -0.0000452 NA -1.7e-06 0 NA NA 10 14300.95 -28581.85 61.40235 4.641190e-14
26 0.1397605 4e-07 -0.0001201 0 -0.0000961 -3.08e-05 0.0000624 NA -4.8e-06 NA 1.64e-05 -5.7e-06 11 14299.18 -28576.32 66.93341 2.921291e-15
5 0.1406261 6e-07 -0.0001493 0 -0.0000090 -6.53e-05 -0.0000857 0e+00 NA 0 NA NA 10 14296.55 -28573.06 70.19982 5.705361e-16
25 0.1402082 8e-07 -0.0001405 0 -0.0000999 -4.30e-05 -0.0000430 -1e-07 NA NA 1.76e-05 -2.6e-06 11 14294.92 -28567.80 75.45446 4.123369e-17
11 0.1399495 7e-07 -0.0001307 0 -0.0000978 -5.69e-05 0.0000099 -1e-07 -2.0e-06 NA 1.73e-05 NA 11 14294.66 -28567.27 75.98128 3.168506e-17
19 0.1400559 5e-07 -0.0001333 0 -0.0000113 -5.77e-05 0.0000174 0e+00 -2.7e-06 NA NA -1.7e-06 11 14293.70 -28565.35 77.90868 1.208721e-17
28 0.1406725 4e-07 -0.0001470 0 -0.0000935 -3.05e-05 -0.0001055 NA NA 0 1.64e-05 -5.4e-06 11 14291.80 -28561.55 81.70210 1.813827e-18
14 0.1402430 4e-07 -0.0001316 0 -0.0000890 -5.86e-05 -0.0000209 NA -3.0e-06 0 1.58e-05 NA 11 14291.46 -28560.87 82.38244 1.290804e-18
22 0.1404322 5e-07 -0.0001410 0 -0.0000103 -5.50e-05 -0.0000416 NA -1.9e-06 0 NA -2.3e-06 11 14290.57 -28559.08 84.17121 5.277578e-19
13 0.1404724 8e-07 -0.0001438 0 -0.0000991 -5.62e-05 -0.0000802 -1e-07 NA 0 1.79e-05 NA 11 14287.28 -28552.52 90.73241 1.984664e-20
21 0.1406404 6e-07 -0.0001497 0 -0.0000094 -5.93e-05 -0.0000874 0e+00 NA 0 NA -1.2e-06 11 14286.20 -28550.36 92.89788 6.721422e-21
7 0.1404979 6e-07 -0.0001445 0 -0.0000092 -6.58e-05 -0.0000590 0e+00 -1.0e-06 0 NA NA 11 14285.96 -28549.88 93.37861 5.285327e-21
27 0.1398914 7e-07 -0.0001279 0 -0.0001008 -3.94e-05 0.0000252 -1e-07 -2.8e-06 NA 1.76e-05 -3.6e-06 12 14284.39 -28544.73 98.52951 4.023185e-22
30 0.1401987 4e-07 -0.0001299 0 -0.0000970 -2.90e-05 -0.0000096 NA -3.5e-06 0 1.69e-05 -6.0e-06 12 14281.20 -28538.34 104.91315 1.653402e-23
29 0.1405076 8e-07 -0.0001446 0 -0.0001019 -4.01e-05 -0.0000847 -1e-07 NA 0 1.82e-05 -3.3e-06 12 14276.97 -28529.89 113.36490 2.416050e-25
15 0.1404182 8e-07 -0.0001418 0 -0.0000991 -5.65e-05 -0.0000689 -1e-07 -4.0e-07 0 1.79e-05 NA 12 14276.70 -28529.34 113.91919 1.831227e-25
23 0.1404716 5e-07 -0.0001432 0 -0.0000099 -5.75e-05 -0.0000518 0e+00 -1.4e-06 0 NA -1.7e-06 12 14275.66 -28527.26 115.99255 6.494083e-26
31 0.1403601 7e-07 -0.0001390 0 -0.0001023 -3.86e-05 -0.0000535 -1e-07 -1.2e-06 0 1.82e-05 -3.7e-06 13 14266.43 -28506.79 136.46456 2.328506e-30

From the dredged model list we can see that the top model is one with no interactions. Run this model on the MCC and 100 trees. Noting that each model uses a unique set of values for DR/ND and a unique tree in the correlation structure.

#In both cases the top model is 1/2/3/4/5/6 no interaction terms. With no models within delta < 4: 

#Run model for DR
MCC.DR.top <- gls(MCC.DR ~ SDi
                         + log(range.size.m2)
                         + bioclim4 #Seasonal variation
                         + residuals.PC1 #Spatial variation
                         + PC1.LIG #Long-term climate variation
                         + NPP,
                correlation = corPagel(0.985, phy = pruned.MCC.tree, fixed = TRUE), 
                data = restricted.data, 
                method = "REML")

saveRDS(MCC.DR.top, 'data/MCC.DR.top.rds')

#Run model for ND
MCC.ND.top <- gls(MCC.ND ~ SDi
                      + log(range.size.m2)
                      + bioclim4 #Seasonal variation
                      + residuals.PC1 #Spatial variation
                      + PC1.LIG #Long-term climate variation
                      + NPP,
                correlation = corPagel(0.9996, phy = pruned.MCC.tree, fixed = TRUE), 
                data = restricted.data, 
                method = "REML")

saveRDS(MCC.ND.top, 'data/MCC.ND.top.rds')
#Run the 100 models for DR and ND using the best model:

#Take the restricted data and make it simpler with just responses and predictors.Note that we join the es.values for the 100 trees
DR.model.data <- lapply(es.list, function(x) { #es.list is a list of ES values calculated earlier
  left_join(restricted.data %>% dplyr::select(binomial, TipLabel, SDi, range.size.m2, bioclim4, residuals.PC1, PC1.LIG, NPP), 
            x %>% as.data.frame() %>% tibble::rownames_to_column() %>% `colnames<-`(c("TipLabel", "ES")), 
            by = "TipLabel")
})

#PGLS needs tiplabel as rowname
DR.model.data <- lapply(DR.model.data, function(x) {
  tibble::column_to_rownames(x, "TipLabel")})

#Prune the trees
pruned.trees<-lapply(passerine.trees, function(x) {
  drop.tip(x,x$tip.label[-match(restricted.data$TipLabel, x$tip.label)])
})

#Use mapply to create a list of PGLS global models
DR.pgls.models <- mcmapply(function(x,y) {
  gls(ES ~ SDi 
         + log(range.size.m2)
         + bioclim4 #Seasonal variation
         + residuals.PC1 #Spatial variation
         + PC1.LIG #Long-term climate variation
         + NPP,
    corPagel(0.985, phy = y, fixed = TRUE), 
    data = x, 
    method = "REML")
}, x = DR.model.data, y = pruned.trees,
SIMPLIFY = FALSE, #Prevents the catastrophic loss of structure, names and attributes
mc.cores = 8) #Specify core number 

saveRDS(DR.pgls.models, "data/DR.pgls.models.rds")

#Now for Node Density:
ND.model.data <- lapply(nd.list, function(x) {
  left_join(restricted.data %>% dplyr::select(binomial, TipLabel, SDi, range.size.m2, bioclim4, residuals.PC1, PC1.LIG, NPP), 
            x %>% as.data.frame() %>% tibble::rownames_to_column() %>% `colnames<-`(c("TipLabel", "ND")), 
            by = "TipLabel")
})

#PGLS needs tiplabel as rowname
ND.model.data <- lapply(ND.model.data, function(x) {
  tibble::column_to_rownames(x, "TipLabel")})

#Use mapply to create a list of PGLS global models
ND.pgls.models <- mcmapply(function(x,y) {
gls(ND ~ SDi 
         + log(range.size.m2)
         + bioclim4 #Seasonal variation
         + residuals.PC1 #Spatial variation
         + PC1.LIG #Long-term climate variation
         + NPP,
    corPagel(0.9996, phy = y, fixed = TRUE), 
    data = x, 
    method = "REML")
}, x = ND.model.data, y = pruned.trees, 
SIMPLIFY = FALSE, #Prevents the catastrophic loss of structure, names and attributes
mc.cores = 8) #Specify core number 

saveRDS(ND.pgls.models, "data/ND.pgls.models.rds")

From these models we can get a distribution of the estimates for the models alongside the 95 % CIs of the MCC model, this will show the variation beween trees, and the variation associated with the top MCC model.

DR.pgls.models <- readRDS("data/DR.pgls.models.rds")
MCC.DR.top <- readRDS('data/MCC.DR.top.rds')

DR.pgls.summary <- lapply(DR.pgls.models, function(x) {
  data.frame(x$coefficients, confint(x)) %>% tibble::rownames_to_column()
})

MCC.DR.summary <- data.frame(MCC.DR.top$coefficients, confint(MCC.DR.top)) %>% tibble::rownames_to_column()

DR.pgls.summary <- bind_rows(DR.pgls.summary) %>% `colnames<-`(c("Parameter", "Estimate", "LCI", "UCI"))

colnames(MCC.DR.summary) <- c("Parameter", "Estimate", "LCI", "UCI")

parameter_names <- c(
                    `bioclim4` = "Temperature Seasonality",
                    `log(range.size.m2)` = "Range Size (log-transformed)",
                    `NPP` = "NPP",
                    `PC1.LIG` = "Long-term Temperature Variation",
                    `residuals.PC1` = "Spatial Temperature Variation",
                    `SDi` = "Sexual Dichromatism"
                    )

DR.plot <-DR.pgls.summary %>% filter(Parameter != "(Intercept)") %>% 
  ggplot(aes(x = Estimate, y = Parameter, fill = Parameter)) +
  geom_density(aes(y = ..scaled..),
               position = position_nudge(y= 0.6))+
  geom_jitter(shape = 21, alpha = 0.5, size = 0.75)+
  geom_vline(xintercept = 0)+
  theme_minimal()+
  theme(axis.text.y = element_blank(),
        legend.position = "none")

DR.plot <- DR.plot + geom_errorbarh(data = MCC.DR.summary %>% filter(Parameter != "(Intercept)"), 
                 aes(xmin = LCI,
                     xmax = UCI, y = Parameter,
                     color = Parameter), 
                 height = 0, show.legend = F, 
                 position = position_nudge(y= -0.75))+
  
  geom_point(data = MCC.DR.summary %>% filter(Parameter != "(Intercept)"), 
             aes(x = Estimate, y = Parameter, fill = Parameter), 
             shape=21, color = "grey20",
             size = 3,
             position = position_nudge(y= -0.75)) +
  scale_y_discrete(expand = c(0.5,.5))+
  facet_wrap(~ Parameter, scales = "free", nrow = 6, labeller = as_labeller(parameter_names))+
  scale_fill_brewer(palette = "Dark2")+
  scale_color_brewer(palette = "Dark2")+
  theme(panel.grid.major.y = element_blank(),
        axis.title.y = element_blank(),
        axis.text.x = element_text(size = 7),
        plot.margin=unit(c(.5,.5,.5,.5),"cm"))+
  ggtitle("DR Models")




#Now plot this for Node Density (ND):  

ND.pgls.models <- readRDS("data/ND.pgls.models.rds")
MCC.ND.top <- readRDS('data/MCC.ND.top.rds')

ND.pgls.summary <- lapply(ND.pgls.models, function(x) {
  data.frame(x$coefficients, confint(x)) %>% tibble::rownames_to_column()
})

MCC.ND.summary <- data.frame(MCC.ND.top$coefficients, confint(MCC.ND.top)) %>% tibble::rownames_to_column()

ND.pgls.summary <- bind_rows(ND.pgls.summary) %>% `colnames<-`(c("Parameter", "Estimate", "LCI", "UCI"))

colnames(MCC.ND.summary) <- c("Parameter", "Estimate", "LCI", "UCI")

ND.plot <-ND.pgls.summary %>% filter(Parameter != "(Intercept)") %>% 
  ggplot(aes(x = Estimate, y = Parameter, fill = Parameter)) +
  geom_density(aes(y = ..scaled..),
               position = position_nudge(y= 0.6))+
  geom_jitter(shape = 21, alpha = 0.5, size = 0.75)+
  geom_vline(xintercept = 0)+
  theme_minimal()+
  theme(axis.text.y = element_blank(),
        legend.position = "none")

ND.plot <- ND.plot + geom_errorbarh(data = MCC.ND.summary %>% filter(Parameter != "(Intercept)"), 
                 aes(xmin = LCI,
                     xmax = UCI, y = Parameter,
                     color = Parameter), 
                 height = 0, show.legend = F, 
                 position = position_nudge(y= -0.75))+
  
  geom_point(data = MCC.ND.summary %>% filter(Parameter != "(Intercept)"), 
             aes(x = Estimate, y = Parameter, fill = Parameter), 
             shape=21, color = "grey20",
             size = 3,
             position = position_nudge(y= -0.75)) +
  scale_y_discrete(expand = c(0.5,0.5))+
  facet_wrap(~ Parameter, scales = "free", nrow = 6, labeller = as_labeller(parameter_names))+
  scale_fill_brewer(palette = "Dark2")+
  scale_color_brewer(palette = "Dark2")+
  theme(panel.grid.major.y = element_blank(),
        axis.title.y = element_blank(),
        axis.text.x = element_text(size = 7),
        plot.margin=unit(c(.5,.5,.5,.5),"cm"))+
  ggtitle("ND Models")

PGLS Models on BAMM estimates

Using the previously generated BAMM dataframe we can now undertake model selection for speciation and extinction. The process is similar to that used on DR and ND, except given the use of a Bayesian approach in BAMM we can make use of varying levels of uncertainty between tips (species) by constructing a weighted model, where the weight is the inverse of the variance, meaning more precise estimates of speciation or extinction at a given tip (species) holds higher weight in the model.

Based on preliminary findings we found that Pagal’s lambda was = 1 and running corPagel lead to problemd of convergence. Therefore the following models are run assuming Brownian Motion with corBrownian.

MCC.BAMM.df <- readRDS('data/MCC.BAMM.df.rds')

#Create model dataframe for use in models
MCC.BAMM.model.data <- left_join(restricted.data %>% 
                                   dplyr::select(binomial, TipLabel, SDi, range.size.m2, bioclim4, residuals.PC1, PC1.LIG, NPP), 
                                 MCC.BAMM.df %>% as.data.frame(), 
                                 by = "TipLabel")

saveRDS(MCC.BAMM.model.data, 'data/MCC.BAMM.model.data.rds')
#Prune tree
pruned.MCC.tree <- drop.tip(MCC.passerine,MCC.passerine$tip.label[-match(MCC.BAMM.model.data$TipLabel, MCC.passerine$tip.label)])

#Set rownames to match tree
rownames(MCC.BAMM.model.data) <- MCC.BAMM.model.data$TipLabel

#Run model for DR
MCC.BAMM.lambda <- gls(mean.lambda ~ SDi
                         + log(range.size.m2)
                         + bioclim4 #Seasonal variation
                         + residuals.PC1 #Spatial variation
                         + PC1.LIG #Long-term climate variation
                         + NPP
                         + SDi*log(range.size.m2)
                         + SDi*bioclim4
                         + SDi*residuals.PC1
                         + SDi*PC1.LIG
                         + SDi*NPP,
                weights = ~ sqrt(var.lambda), #sqrt to account for overdispersedskewed variance distribution
                correlation = corBrownian(phy = pruned.MCC.tree), 
                data = MCC.BAMM.model.data, 
                method = "REML")

#Run model for DR
MCC.BAMM.mu <- gls(mean.mu ~ SDi
                         + log(range.size.m2)
                         + bioclim4 #Seasonal variation
                         + residuals.PC1 #Spatial variation
                         + PC1.LIG #Long-term climate variation
                         + NPP
                         + SDi*log(range.size.m2)
                         + SDi*bioclim4
                         + SDi*residuals.PC1
                         + SDi*PC1.LIG
                         + SDi*NPP,
                weights = ~ sqrt(var.mu),
                correlation = corBrownian(phy = pruned.MCC.tree), 
                data = MCC.BAMM.model.data, 
                method = "REML")


#Dredge the global MCC models

#Set up cluster
cores<-8
clusterType <- if(length(find.package("snow", quiet = TRUE))) "SOCK" else "PSOCK"
clust <- try(makeCluster(getOption("cl.cores", cores), type = clusterType))
myclust<-clust

#Export data and packages to cluster
clusterExport(myclust, c("MCC.BAMM.model.data"), envir=environment())
clusterExport(myclust, c("pruned.MCC.tree"), envir=environment())
clusterEvalQ(myclust, library(nlme))
clusterEvalQ(myclust, library(ape))
clusterEvalQ(myclust, library(MuMIn))

#Dredged models:

dredged.MCC.lambda <- pdredge(MCC.BAMM.lambda, fixed = c("SDi", "log(range.size.m2)", "bioclim4", "residuals.PC1", "PC1.LIG", "NPP"), cluster=myclust)

saveRDS(dredged.MCC.lambda, "data/dredged.MCC.lambda.rds")

dredged.MCC.mu <- pdredge(MCC.BAMM.mu, fixed = c("SDi", "log(range.size.m2)", "bioclim4", "residuals.PC1", "PC1.LIG", "NPP"), cluster=myclust)

saveRDS(dredged.MCC.mu, "data/dredged.MCC.mu.rds")

Table S6: The dredged models for BAMM speciation and BAMM extinction both show the top model is one with no interactions , with \(\delta AICc > 4\). This is the same situation as DR/ND models (see above).

dredged.MCC.lambda <- readRDS("data/dredged.MCC.lambda.rds")
dredged.MCC.mu <- readRDS("data/dredged.MCC.mu.rds")
kable(dredged.MCC.lambda, "html", caption = "BAMM-Speciation Dredge Table") %>%
  kable_styling() %>%
  scroll_box(width = "100%", height = "500px")
BAMM-Speciation Dredge Table
(Intercept) bioclim4 log(range.size.m2) NPP PC1.LIG residuals.PC1 SDi bioclim4:SDi log(range.size.m2):SDi NPP:SDi PC1.LIG:SDi residuals.PC1:SDi df logLik AICc delta weight
0 -2.109226 -9.0e-07 -0.0002225 0 -0.0000487 0.0000043 -0.0000143 NA NA NA NA NA 8 6820.148 -13624.27 0.00000 9.998911e-01
16 -2.109176 -9.0e-07 -0.0002237 0 -0.0000508 0.0000683 -0.0000199 NA NA NA NA -1.38e-05 9 6811.014 -13604.00 20.27447 3.957376e-05
8 -2.109265 -1.0e-06 -0.0002207 0 -0.0000969 0.0000114 -0.0000161 NA NA NA 1.04e-05 NA 9 6810.971 -13603.91 20.35993 3.791842e-05
2 -2.109396 -9.0e-07 -0.0002160 0 -0.0000493 0.0000040 0.0000283 NA -1.6e-06 NA NA NA 9 6810.770 -13603.51 20.76242 3.100641e-05
1 -2.109089 -1.2e-06 -0.0002250 0 -0.0000511 0.0000046 -0.0000299 1e-07 NA NA NA NA 9 6806.364 -13594.70 29.57313 3.786430e-07
4 -2.108967 -9.0e-07 -0.0002264 0 -0.0000479 0.0000082 -0.0000514 NA NA 0 NA NA 9 6803.214 -13588.40 35.87346 1.622295e-08
24 -2.109217 -1.0e-06 -0.0002215 0 -0.0001185 0.0000950 -0.0000239 NA NA NA 1.46e-05 -1.74e-05 10 6801.896 -13583.75 40.51732 1.591214e-09
18 -2.109482 -9.0e-07 -0.0002120 0 -0.0000517 0.0000695 0.0000568 NA -2.9e-06 NA NA -1.41e-05 10 6801.644 -13583.25 41.02023 1.237434e-09
10 -2.109495 -1.0e-06 -0.0002119 0 -0.0000984 0.0000111 0.0000414 NA -2.2e-06 NA 1.06e-05 NA 10 6801.596 -13583.15 41.11709 1.178937e-09
17 -2.108913 -1.4e-06 -0.0002286 0 -0.0000555 0.0000827 -0.0000498 1e-07 NA NA NA -1.68e-05 10 6797.294 -13574.55 49.72130 1.596283e-11
9 -2.109192 -1.1e-06 -0.0002221 0 -0.0000944 0.0000110 -0.0000239 0e+00 NA NA 9.60e-06 NA 10 6797.223 -13574.41 49.86312 1.487013e-11
3 -2.109506 -1.3e-06 -0.0002080 0 -0.0000535 0.0000038 0.0000821 1e-07 -4.5e-06 NA NA NA 10 6797.094 -13574.15 50.12135 1.306895e-11
20 -2.108936 -9.0e-07 -0.0002273 0 -0.0000499 0.0000705 -0.0000544 NA NA 0 NA -1.34e-05 10 6794.078 -13568.12 56.15218 6.407108e-13
12 -2.108992 -1.0e-06 -0.0002247 0 -0.0000977 0.0000158 -0.0000554 NA NA 0 1.08e-05 NA 10 6794.041 -13568.04 56.22755 6.170159e-13
6 -2.108998 -9.0e-07 -0.0002253 0 -0.0000480 0.0000081 -0.0000438 NA -3.0e-07 0 NA NA 10 6793.855 -13567.67 56.59936 5.123386e-13
26 -2.109645 -1.0e-06 -0.0002051 0 -0.0001220 0.0000975 0.0000831 NA -4.0e-06 NA 1.50e-05 -1.80e-05 11 6792.534 -13563.02 61.24819 5.012765e-14
5 -2.108812 -1.2e-06 -0.0002292 0 -0.0000504 0.0000086 -0.0000693 1e-07 NA 0 NA NA 10 6789.433 -13558.83 65.44243 6.156164e-15
25 -2.109027 -1.3e-06 -0.0002252 0 -0.0001141 0.0001021 -0.0000445 1e-07 NA NA 1.29e-05 -1.91e-05 11 6788.188 -13554.33 69.94123 6.492436e-16
19 -2.109734 -1.5e-06 -0.0001945 0 -0.0000611 0.0000941 0.0001749 2e-07 -9.1e-06 NA NA -1.95e-05 11 6788.081 -13554.12 70.15475 5.835015e-16
11 -2.109552 -1.2e-06 -0.0002074 0 -0.0000951 0.0000101 0.0000733 1e-07 -3.9e-06 NA 9.30e-06 NA 11 6787.953 -13553.86 70.41057 5.134414e-16
28 -2.108963 -1.0e-06 -0.0002253 0 -0.0001189 0.0000977 -0.0000606 NA NA 0 1.48e-05 -1.71e-05 11 6784.963 -13547.88 76.38995 2.582768e-17
22 -2.109125 -9.0e-07 -0.0002204 0 -0.0000506 0.0000710 -0.0000083 NA -1.7e-06 0 NA -1.37e-05 11 6784.727 -13547.41 76.86170 2.040073e-17
14 -2.109083 -1.0e-06 -0.0002214 0 -0.0000982 0.0000156 -0.0000334 NA -8.0e-07 0 1.09e-05 NA 11 6784.683 -13547.32 76.95041 1.951567e-17
21 -2.108655 -1.4e-06 -0.0002325 0 -0.0000548 0.0000854 -0.0000867 1e-07 NA 0 NA -1.65e-05 11 6780.361 -13538.68 85.59503 2.589570e-19
13 -2.108909 -1.2e-06 -0.0002264 0 -0.0000949 0.0000153 -0.0000645 0e+00 NA 0 9.90e-06 NA 11 6780.294 -13538.54 85.72902 2.421769e-19
7 -2.109119 -1.3e-06 -0.0002173 0 -0.0000521 0.0000078 0.0000103 1e-07 -3.0e-06 0 NA NA 11 6780.178 -13538.31 85.95984 2.157800e-19
27 -2.109821 -1.5e-06 -0.0001922 0 -0.0001185 0.0001127 0.0001733 1e-07 -8.8e-06 NA 1.27e-05 -2.17e-05 12 6778.973 -13533.89 90.37795 2.369350e-20
30 -2.109279 -1.0e-06 -0.0002136 0 -0.0001213 0.0000993 0.0000163 NA -2.8e-06 0 1.51e-05 -1.76e-05 12 6775.618 -13527.18 97.08814 8.270481e-22
29 -2.108760 -1.3e-06 -0.0002292 0 -0.0001143 0.0001051 -0.0000828 1e-07 NA 0 1.31e-05 -1.88e-05 12 6771.256 -13518.46 105.81172 1.054947e-23
23 -2.109417 -1.5e-06 -0.0002023 0 -0.0000598 0.0000946 0.0001146 2e-07 -7.8e-06 0 NA -1.90e-05 12 6771.163 -13518.27 105.99910 9.605952e-24
15 -2.109145 -1.2e-06 -0.0002171 0 -0.0000953 0.0000146 -0.0000028 0e+00 -2.4e-06 0 9.70e-06 NA 12 6771.041 -13518.03 106.24362 8.500498e-24
31 -2.109486 -1.5e-06 -0.0002005 0 -0.0001180 0.0001135 0.0001092 1e-07 -7.4e-06 0 1.29e-05 -2.11e-05 13 6762.057 -13498.05 126.21895 3.907117e-28
kable(dredged.MCC.mu, "html", caption = "BAMM-Extinction Dredge Table") %>%
  kable_styling() %>%
  scroll_box(width = "100%", height = "500px")
BAMM-Extinction Dredge Table
(Intercept) bioclim4 log(range.size.m2) NPP PC1.LIG residuals.PC1 SDi bioclim4:SDi log(range.size.m2):SDi NPP:SDi PC1.LIG:SDi residuals.PC1:SDi df logLik AICc delta weight
0 -4.355741 -1.3e-06 -0.0004220 0e+00 -0.0002276 -0.0000745 0.0000239 NA NA NA NA NA 8 2587.288 -5158.551 0.00000 9.997936e-01
16 -4.355884 -1.2e-06 -0.0004231 0e+00 -0.0002369 0.0000641 0.0000107 NA NA NA NA -2.98e-05 9 2578.877 -5139.723 18.82838 8.154169e-05
8 -4.355680 -1.3e-06 -0.0004247 0e+00 -0.0001838 -0.0000819 0.0000266 NA NA NA -8.50e-06 NA 9 2578.674 -5139.317 19.23423 6.656558e-05
2 -4.356963 -1.3e-06 -0.0003779 0e+00 -0.0002336 -0.0000760 0.0002793 NA -9.60e-06 NA NA NA 9 2578.528 -5139.025 19.52584 5.753445e-05
1 -4.355554 -1.7e-06 -0.0004251 0e+00 -0.0002285 -0.0000749 0.0000004 1e-07 NA NA NA NA 9 2574.173 -5130.314 28.23710 7.384157e-07
4 -4.354974 -1.3e-06 -0.0004343 -1e-07 -0.0002234 -0.0000731 -0.0000587 NA NA 0 NA NA 9 2570.982 -5123.933 34.61801 3.038819e-08
24 -4.355849 -1.2e-06 -0.0004245 0e+00 -0.0002139 0.0000564 0.0000125 NA NA NA -4.40e-06 -2.90e-05 10 2570.271 -5120.505 38.04631 5.473417e-09
18 -4.357443 -1.1e-06 -0.0003672 0e+00 -0.0002451 0.0000706 0.0003338 NA -1.21e-05 NA NA -3.16e-05 10 2570.130 -5120.223 38.32863 4.752851e-09
10 -4.356868 -1.2e-06 -0.0003818 0e+00 -0.0001920 -0.0000829 0.0002739 NA -9.30e-06 NA -8.10e-06 NA 10 2569.914 -5119.791 38.76025 3.830269e-09
17 -4.355503 -2.0e-06 -0.0004301 0e+00 -0.0002406 0.0000909 -0.0000433 2e-07 NA NA NA -3.58e-05 10 2565.829 -5111.620 46.93151 6.439636e-11
9 -4.355419 -1.8e-06 -0.0004295 0e+00 -0.0001730 -0.0000844 -0.0000034 1e-07 NA NA -1.08e-05 NA 10 2565.588 -5111.139 47.41245 5.063225e-11
3 -4.357282 -2.0e-06 -0.0003593 0e+00 -0.0002385 -0.0000775 0.0003758 2e-07 -1.49e-05 NA NA NA 10 2565.518 -5110.999 47.55279 4.720112e-11
20 -4.355085 -1.1e-06 -0.0004360 -1e-07 -0.0002327 0.0000686 -0.0000759 NA NA 0 NA -3.05e-05 10 2562.574 -5105.111 53.44028 2.485999e-12
12 -4.354937 -1.2e-06 -0.0004364 -1e-07 -0.0001843 -0.0000798 -0.0000544 NA NA 0 -7.60e-06 NA 10 2562.368 -5104.698 53.85291 2.022555e-12
6 -4.355914 -1.2e-06 -0.0004017 -1e-07 -0.0002281 -0.0000743 0.0001311 NA -6.80e-06 0 NA NA 10 2562.236 -5104.434 54.11718 1.772201e-12
26 -4.357390 -1.1e-06 -0.0003691 0e+00 -0.0002265 0.0000643 0.0003302 NA -1.20e-05 NA -3.60e-06 -3.09e-05 11 2561.526 -5101.007 57.54451 3.193580e-13
5 -4.354737 -1.7e-06 -0.0004382 -1e-07 -0.0002242 -0.0000736 -0.0000879 1e-07 NA 0 NA NA 10 2557.870 -5095.703 62.84860 2.251695e-14
19 -4.358289 -2.6e-06 -0.0003247 0e+00 -0.0002597 0.0001253 0.0005543 3e-07 -2.41e-05 NA NA -4.41e-05 11 2557.246 -5092.447 66.10399 4.421925e-15
25 -4.355407 -2.1e-06 -0.0004331 0e+00 -0.0002004 0.0000796 -0.0000449 2e-07 NA NA -7.80e-06 -3.48e-05 11 2557.245 -5092.444 66.10718 4.414887e-15
11 -4.357205 -2.1e-06 -0.0003614 0e+00 -0.0001787 -0.0000879 0.0003868 2e-07 -1.55e-05 NA -1.18e-05 NA 11 2556.938 -5091.830 66.72172 3.246924e-15
28 -4.355066 -1.1e-06 -0.0004369 -1e-07 -0.0002154 0.0000628 -0.0000737 NA NA 0 -3.30e-06 -2.98e-05 11 2553.970 -5085.894 72.65695 1.669761e-16
22 -4.356379 -1.1e-06 -0.0003913 -1e-07 -0.0002395 0.0000731 0.0001834 NA -9.40e-06 0 NA -3.18e-05 11 2553.839 -5085.633 72.91876 1.464880e-16
14 -4.355849 -1.2e-06 -0.0004048 -1e-07 -0.0001901 -0.0000806 0.0001292 NA -6.60e-06 0 -7.40e-06 NA 11 2553.622 -5085.199 73.35222 1.179447e-16
21 -4.354615 -2.1e-06 -0.0004445 -1e-07 -0.0002364 0.0000981 -0.0001408 2e-07 NA 0 NA -3.70e-05 11 2549.534 -5077.022 81.52954 1.976951e-18
13 -4.354629 -1.8e-06 -0.0004420 -1e-07 -0.0001726 -0.0000824 -0.0000896 1e-07 NA 0 -1.01e-05 NA 11 2549.286 -5076.526 82.02573 1.542585e-18
7 -4.356245 -2.0e-06 -0.0003832 -1e-07 -0.0002330 -0.0000758 0.0002278 2e-07 -1.21e-05 0 NA NA 11 2549.225 -5076.405 82.14605 1.452519e-18
27 -4.358212 -2.7e-06 -0.0003270 0e+00 -0.0002158 0.0001133 0.0005585 4e-07 -2.43e-05 NA -8.50e-06 -4.31e-05 12 2548.664 -5073.274 85.27712 3.035414e-19
30 -4.356346 -1.1e-06 -0.0003927 -1e-07 -0.0002249 0.0000682 0.0001817 NA -9.30e-06 0 -2.80e-06 -3.12e-05 12 2545.236 -5066.417 92.13388 9.846754e-21
23 -4.357251 -2.6e-06 -0.0003486 -1e-07 -0.0002541 0.0001270 0.0004062 3e-07 -2.13e-05 0 NA -4.41e-05 12 2540.954 -5057.854 100.69765 1.360485e-22
29 -4.354544 -2.1e-06 -0.0004469 -1e-07 -0.0002010 0.0000881 -0.0001407 2e-07 NA 0 -6.90e-06 -3.61e-05 12 2540.950 -5057.845 100.70611 1.354739e-22
15 -4.356211 -2.1e-06 -0.0003843 -1e-07 -0.0001773 -0.0000856 0.0002434 2e-07 -1.27e-05 0 -1.10e-05 NA 12 2540.644 -5057.235 101.31622 9.985545e-23
31 -4.357207 -2.6e-06 -0.0003500 -1e-07 -0.0002146 0.0001161 0.0004136 4e-07 -2.16e-05 0 -7.70e-06 -4.32e-05 13 2532.372 -5038.680 119.87105 9.337735e-27

Given that we are running models weighted according to the variance of the response (\(\lambda_{BAMM}\) and \(\mu_{BAMM}\)), we can check to see whether a weighted model isw favourable to an unweighted model using an anova to compare AIC values.

#Run the top model for the MCC
#Run model for ND
MCC.Lambda.top <- gls(mean.lambda ~ SDi
                      + log(range.size.m2)
                      + bioclim4 #Seasonal variation
                      + residuals.PC1 #Spatial variation
                      + PC1.LIG #Long-term climate variation
                      + NPP,
                weights = ~ sqrt(var.lambda),
                correlation = corBrownian(phy = pruned.MCC.tree), 
                data = MCC.BAMM.model.data, 
                method = "REML")
saveRDS(MCC.Lambda.top, 'data/MCC.Lambda.top.rds')

MCC.Mu.top <- gls(mean.mu ~ SDi
                      + log(range.size.m2)
                      + bioclim4 #Seasonal variation
                      + residuals.PC1 #Spatial variation
                      + PC1.LIG #Long-term climate variation
                      + NPP,
                weights = ~ sqrt(var.mu),
                correlation = corBrownian(phy = pruned.MCC.tree), 
                data = MCC.BAMM.model.data, 
                method = "REML")

saveRDS(MCC.Mu.top, 'data/MCC.Mu.top.rds')


#We can also see how the models look without the weightings: 
MCC.Lambda.top.unweighted <- gls(mean.lambda ~ SDi
                      + log(range.size.m2)
                      + bioclim4 #Seasonal variation
                      + residuals.PC1 #Spatial variation
                      + PC1.LIG #Long-term climate variation
                      + NPP,
                correlation = corBrownian(phy = pruned.MCC.tree), 
                data = MCC.BAMM.model.data, 
                method = "REML")
saveRDS(MCC.Lambda.top.unweighted, 'data/MCC.Lambda.top.unweighted.rds')


MCC.Mu.top.unweighted <- gls(mean.mu ~ SDi
                      + log(range.size.m2)
                      + bioclim4 #Seasonal variation
                      + residuals.PC1 #Spatial variation
                      + PC1.LIG #Long-term climate variation
                      + NPP,
                correlation = corBrownian(phy = pruned.MCC.tree), 
                data = MCC.BAMM.model.data, 
                method = "REML")
saveRDS(MCC.Mu.top.unweighted, 'data/MCC.Mu.top.unweighted.rds')

Table S6: Models with a set of weights for the response based on the inverse of the variance of the posterior distribution in \(\lambda_{BAMM}\) and \(\mu_{BAMM}\) have much lower AIC values, indicating that the weighting scheme improves the fit of the model. In any case, we found that an unweighted model did not change the main conclusions drawn in that they did not generate any statistically significant relationships.

MCC.Lambda.top <- readRDS('data/MCC.Lambda.top.rds')
MCC.Lambda.top.unweighted <- readRDS('data/MCC.Lambda.top.unweighted.rds')
MCC.Mu.top <- readRDS('data/MCC.Mu.top.rds')
MCC.Mu.top.unweighted <- readRDS('data/MCC.Mu.top.unweighted.rds')

anova(MCC.Lambda.top, MCC.Lambda.top.unweighted) %>% pander(split.table = Inf)
  call Model df AIC BIC logLik
MCC.Lambda.top gls(model = mean.lambda ~ SDi + log(range.size.m2) + bioclim4 + residuals.PC1 + PC1.LIG + NPP, data = MCC.BAMM.model.data, correlation = corBrownian(phy = pruned.MCC.tree), weights = ~sqrt(var.lambda), method = “REML”) 1 8 -13624 -13571 6820
MCC.Lambda.top.unweighted gls(model = mean.lambda ~ SDi + log(range.size.m2) + bioclim4 + residuals.PC1 + PC1.LIG + NPP, data = MCC.BAMM.model.data, correlation = corBrownian(phy = pruned.MCC.tree), method = “REML”) 2 8 -11788 -11735 5902
anova(MCC.Mu.top, MCC.Mu.top.unweighted) %>% pander(split.table = Inf)
  call Model df AIC BIC logLik
MCC.Mu.top gls(model = mean.mu ~ SDi + log(range.size.m2) + bioclim4 + residuals.PC1 + PC1.LIG + NPP, data = MCC.BAMM.model.data, correlation = corBrownian(phy = pruned.MCC.tree), weights = ~sqrt(var.mu), method = “REML”) 1 8 -5159 -5105 2587
MCC.Mu.top.unweighted gls(model = mean.mu ~ SDi + log(range.size.m2) + bioclim4 + residuals.PC1 + PC1.LIG + NPP, data = MCC.BAMM.model.data, correlation = corBrownian(phy = pruned.MCC.tree), method = “REML”) 2 8 -5320 -5266 2668

We can run the top (no interactions) model on the 100 trees, each with unique estimates of speciation and extinction from the BAMM runs.

#Read in the BAMM data for the 100 trees
BAMM.df <- readRDS('data/BAMM.df.rds')

#Take the restricted data and make it simpler with just responses and predictors.Note that we join the BAMM for the 100 trees

BAMM.model.data <- lapply(BAMM.df, function(x) { #es.list is a list of ES values calculated earlier
  left_join(restricted.data %>% dplyr::select(binomial, TipLabel, SDi, range.size.m2, bioclim4, residuals.PC1, PC1.LIG, NPP), 
            x %>% as.data.frame(), 
            by = "TipLabel")
})

#PGLS needs tiplabel as rowname
BAMM.model.data <- lapply(BAMM.model.data, function(x) {
  tibble::column_to_rownames(x, "TipLabel")})

#Prune the trees
pruned.trees<-lapply(passerine.trees, function(x) {
  drop.tip(x,x$tip.label[-match(restricted.data$TipLabel, x$tip.label)])
})

#Use mapply to create a list of PGLS global models
BAMM.lambda.pgls.models <- mcmapply(function(x,y) {
  gls(mean.lambda ~ SDi 
         + log(range.size.m2)
         + bioclim4 #Seasonal variation
         + residuals.PC1 #Spatial variation
         + PC1.LIG #Long-term climate variation
         + NPP,
    weights = ~ sqrt(var.lambda),      
    corBrownian(phy = y), 
    data = x, 
    method = "REML")
}, x = BAMM.model.data, y = pruned.trees,
SIMPLIFY = FALSE, #Prevents the catastrophic loss of structure, names and attributes
mc.cores = 8) #Specify core number 

saveRDS(BAMM.lambda.pgls.models, "data/BAMM.lambda.pgls.models.rds")

#Use mapply to create a list of PGLS global models
BAMM.mu.pgls.models <- mcmapply(function(x,y) {
  gls(mean.mu ~ SDi 
         + log(range.size.m2)
         + bioclim4 #Seasonal variation
         + residuals.PC1 #Spatial variation
         + PC1.LIG #Long-term climate variation
         + NPP,
    weights = ~ sqrt(var.mu),      
    corBrownian(phy = y), 
    data = x, 
    method = "REML")
}, x = BAMM.model.data, y = pruned.trees,
SIMPLIFY = FALSE, #Prevents the catastrophic loss of structure, names and attributes
mc.cores = 8) #Specify core number 

saveRDS(BAMM.mu.pgls.models, "data/BAMM.mu.pgls.models.rds")
BAMM.lambda.pgls.models <- readRDS("data/BAMM.lambda.pgls.models.rds")
MCC.Lambda.top <- readRDS('data/MCC.Lambda.top.rds')

BAMM.lambda.pgls.summary <- lapply(BAMM.lambda.pgls.models, function(x) {
  data.frame(x$coefficients, confint(x)) %>% tibble::rownames_to_column()
})

parameter_names <- c(
                    `bioclim4` = "Temperature Seasonality",
                    `log(range.size.m2)` = "Range Size (log-transformed)",
                    `NPP` = "NPP",
                    `PC1.LIG` = "Long-term Temperature Variation",
                    `residuals.PC1` = "Spatial Temperature Variation",
                    `SDi` = "Sexual Dichromatism"
                    )

MCC.lambda.summary <- data.frame(MCC.Lambda.top$coefficients, confint(MCC.Lambda.top)) %>% tibble::rownames_to_column()

BAMM.lambda.pgls.summary <- bind_rows(BAMM.lambda.pgls.summary) %>% `colnames<-`(c("Parameter", "Estimate", "LCI", "UCI"))

colnames(MCC.lambda.summary) <- c("Parameter", "Estimate", "LCI", "UCI")

  # for (x in unique(BAMM.lambda.pgls.summary$Parameter)[2:7]){
  #   filter(Parameter == x & between(Estimate, left = as.numeric(hpd.Lambda.top[1,x]), right = as.numeric(hpd.Lambda.top[2,x])))
  #   } 


remove_outliers <- function(x, na.rm = TRUE, ...) {
  qnt <- quantile(x, probs=c(.25, .75), na.rm = na.rm, ...)
  H <- 1.5 * IQR(x, na.rm = na.rm)
  y <- x
  y[x < (qnt[1] - H)] <- NA
  y[x > (qnt[2] + H)] <- NA
  y
}

BAMM.lambda.pgls.summary.RO <- dcast(BAMM.lambda.pgls.summary %>% filter(Parameter != "(Intercept)"), Estimate ~ Parameter, value.var = "Estimate")
BAMM.lambda.pgls.summary.RO$Estimate <- NULL
BAMM.lambda.pgls.summary.RO <- sapply(BAMM.lambda.pgls.summary.RO, function(x) {
  remove_outliers(x, na.rm = T)})
BAMM.lambda.pgls.summary.RO <-melt(BAMM.lambda.pgls.summary.RO) %>% na.omit()
BAMM.lambda.pgls.summary.RO$Var1 <- NULL
colnames(BAMM.lambda.pgls.summary.RO) <- c("Parameter", "Estimate")

BAMM.lambda.plot <- BAMM.lambda.pgls.summary.RO %>%
  ggplot(aes(x = Estimate, y = Parameter, fill = Parameter)) +
  geom_density(aes(y = ..scaled..),
               position = position_nudge(y= 0.6))+
  geom_jitter(shape = 21, alpha = 0.5, size = 0.75)+
  geom_vline(xintercept = 0)+
  theme_minimal()+
  theme(axis.text.y = element_blank(),
        legend.position = "none")

BAMM.lambda.plot <- BAMM.lambda.plot + geom_errorbarh(data = MCC.lambda.summary %>% filter(Parameter != "(Intercept)"), 
                 aes(xmin = LCI,
                     xmax = UCI, y = Parameter,
                     color = Parameter), 
                 height = 0, show.legend = F, 
                 position = position_nudge(y= -0.75))+
  
  geom_point(data = MCC.lambda.summary %>% filter(Parameter != "(Intercept)"), 
             aes(x = Estimate, y = Parameter, fill = Parameter), 
             shape=21, color = "grey20",
             size = 3,
             position = position_nudge(y= -0.75)) +
  scale_y_discrete(expand = c(0.5,0.5))+
  facet_wrap(~ Parameter, scales = "free", nrow = 7, labeller = as_labeller(parameter_names))+
  scale_fill_brewer(palette = "Dark2")+
  scale_color_brewer(palette = "Dark2")+
  theme(panel.grid.major.y = element_blank(),
        axis.title.y = element_blank(),
        axis.text.x = element_text(size = 7),
        plot.margin=unit(c(.5,.5,.5,.5),"cm"))+
  ggtitle("BAMM Speciation")


#Figure for extinction

BAMM.mu.pgls.models <- readRDS("data/BAMM.mu.pgls.models.rds")
MCC.Mu.top <- readRDS('data/MCC.Mu.top.rds')

BAMM.mu.pgls.summary <- lapply(BAMM.mu.pgls.models, function(x) {
  data.frame(x$coefficients, confint(x)) %>% tibble::rownames_to_column()
})

MCC.mu.summary <- data.frame(MCC.Mu.top$coefficients, confint(MCC.Mu.top)) %>% tibble::rownames_to_column()

BAMM.mu.pgls.summary <- bind_rows(BAMM.mu.pgls.summary) %>% `colnames<-`(c("Parameter", "Estimate", "LCI", "UCI"))

colnames(MCC.mu.summary) <- c("Parameter", "Estimate", "LCI", "UCI")

#Get rid of outliers
BAMM.mu.pgls.summary.RO <- dcast(BAMM.mu.pgls.summary %>% select(Parameter,Estimate) %>% filter(Parameter != "(Intercept)"), Estimate ~ Parameter, value.var = "Estimate")
BAMM.mu.pgls.summary.RO$Estimate <- NULL
BAMM.mu.pgls.summary.RO <- sapply(BAMM.mu.pgls.summary.RO, function(x) {
  remove_outliers(x, na.rm = T)})
BAMM.mu.pgls.summary.RO <-melt(BAMM.mu.pgls.summary.RO) %>% na.omit()
BAMM.mu.pgls.summary.RO$Var1 <- NULL
colnames(BAMM.mu.pgls.summary.RO) <- c("Parameter", "Estimate")

BAMM.mu.plot <-BAMM.mu.pgls.summary.RO %>%
  ggplot(aes(x = Estimate, y = Parameter, fill = Parameter)) +
  geom_density(aes(y = ..scaled..),
               position = position_nudge(y= 0.6))+
  geom_jitter(shape = 21, alpha = 0.5, size = 0.75)+
  geom_vline(xintercept = 0)+
  theme_minimal()+
  theme(axis.text.y = element_blank(),
        legend.position = "none")

BAMM.mu.plot <- BAMM.mu.plot + geom_errorbarh(data = MCC.mu.summary %>% filter(Parameter != "(Intercept)"), 
                 aes(xmin = LCI,
                     xmax = UCI, y = Parameter,
                     color = Parameter), 
                 height = 0, show.legend = F, 
                 position = position_nudge(y= -0.75))+
  
  geom_point(data = MCC.mu.summary %>% filter(Parameter != "(Intercept)"), 
             aes(x = Estimate, y = Parameter, fill = Parameter), 
             shape=21, color = "grey20",
             size = 3,
             position = position_nudge(y= -0.75)) +
  scale_y_discrete(expand = c(0.5,0.5))+
  facet_wrap(~ Parameter, scales = "free", nrow = 7, labeller = as_labeller(parameter_names))+
  scale_fill_brewer(palette = "Dark2")+
  scale_color_brewer(palette = "Dark2")+
  theme(panel.grid.major.y = element_blank(),
        axis.title.y = element_blank(),
        axis.text.x = element_text(size = 7),
        plot.margin=unit(c(.5,.5,0.5,.5),"cm"))+
  ggtitle("BAMM Extinction")

PGLS Model Summary

grid.arrange(
  symmetrise_scale(DR.plot, "x"),
  symmetrise_scale(ND.plot, "x"),
  symmetrise_scale(BAMM.lambda.plot, "x"),
  symmetrise_scale(BAMM.mu.plot, "x"),
  nrow = 1
)

Figure S8: This figure is the same basic figure as seen in the manuscript (FIgure 1), it provides model estimates for four response variables across 100 random trees alongside the MCC tree.

Table S7: The estimates of the MCC models plotted above are based on the following data tables.

MCC.DR.summary %>% filter(Parameter != "(Intercept)") %>% pander(caption = "MCC DR Estimates")
MCC DR Estimates
Parameter Estimate LCI UCI
SDi -0.001279 -0.003324 0.0007665
log(range.size.m2) -0.006578 -0.01078 -0.002381
bioclim4 2.141e-05 -3.369e-05 7.65e-05
residuals.PC1 -0.00243 -0.007562 0.002703
PC1.LIG 0.001586 -0.003676 0.006849
NPP 2.937e-07 -1.321e-06 1.909e-06
MCC.ND.summary %>% filter(Parameter != "(Intercept)") %>% pander(caption = "MCC ND Estimates")
MCC ND Estimates
Parameter Estimate LCI UCI
SDi -5.745e-05 -0.0001239 8.957e-06
log(range.size.m2) -0.0001462 -0.0002827 -9.773e-06
bioclim4 4.523e-07 -1.426e-06 2.33e-06
residuals.PC1 -6.469e-05 -0.0002368 0.0001074
PC1.LIG -1.017e-05 -0.0001763 0.000156
NPP -9.069e-09 -6.454e-08 4.64e-08
MCC.lambda.summary %>% filter(Parameter != "(Intercept)") %>% pander(caption = "MCC BAMM Speciation Estimates")
MCC BAMM Speciation Estimates
Parameter Estimate LCI UCI
SDi -1.43e-05 -0.0002574 0.0002288
log(range.size.m2) -0.0002225 -0.00069 0.000245
bioclim4 -9.418e-07 -7.242e-06 5.359e-06
residuals.PC1 4.311e-06 -0.0005798 0.0005884
PC1.LIG -4.874e-05 -0.0005747 0.0004772
NPP -2.09e-08 -2.282e-07 1.864e-07
MCC.mu.summary %>% filter(Parameter != "(Intercept)") %>% pander(caption = "MCC BAMM Extinction Estimates")
MCC BAMM Extinction Estimates
Parameter Estimate LCI UCI
SDi 2.39e-05 -0.0004098 0.0004576
log(range.size.m2) -0.000422 -0.00133 0.0004858
bioclim4 -1.289e-06 -1.401e-05 1.143e-05
residuals.PC1 -7.453e-05 -0.001231 0.001082
PC1.LIG -0.0002276 -0.001284 0.0008287
NPP -2.024e-08 -4.178e-07 3.773e-07

From the estimates of 100 trees we can calculate HPD intervals. Note that these do no account for the CI of each estimate. Thus it is a more a representation of tree uncertainty than model uncertainty.

Table S8: Highest posterior density intervals are calculated for the above figure (ie; models using RGB values of sexual dichromatism). The HPD range is determined using the hdi function with a 95 % credible interval. These intervals do not take into account the variance associated with each interval and thus are not an estimate of model precision. Intervals not overlapping zero suggest that 95 % of trees from the posterior generate a model estimate for the given parameter that in the same direction (+ or -).

hpd.DR.top <- list()
for(x in unique(DR.pgls.summary$Parameter)) {
hpd.DR.top[[x]] = hdi(DR.pgls.summary %>% filter(Parameter == x) %>% dplyr::select("Estimate"))
}
hpd.DR.top <- bind_rows(hpd.DR.top) %>% `rownames<-`(c("Lower", "Upper")) %>% dplyr::select(-"(Intercept)")


saveRDS(hpd.DR.top, 'data/hpd.DR.top.rds')

#For ND
hpd.ND.top <- list()
for(x in unique(ND.pgls.summary$Parameter)) {
hpd.ND.top[[x]] = hdi(ND.pgls.summary %>% filter(Parameter == x) %>% dplyr::select(Estimate))
}
hpd.ND.top <- bind_rows(hpd.ND.top) %>% `rownames<-`(c("Lower", "Upper")) %>% dplyr::select(-"(Intercept)") 

saveRDS(hpd.ND.top, 'data/hpd.ND.top.rds')

hpd.Lambda.top <- list()
for(x in unique(BAMM.lambda.pgls.summary$Parameter)) {
hpd.Lambda.top[[x]] = hdi(BAMM.lambda.pgls.summary %>% filter(Parameter == x) %>% dplyr::select(Estimate))
}
hpd.Lambda.top <- bind_rows(hpd.Lambda.top) %>% `rownames<-`(c("Lower", "Upper")) %>% dplyr::select(-"(Intercept)") 


saveRDS(hpd.Lambda.top, 'data/hpd.Lambda.top.rds')

hpd.Mu.top <- list()
for(x in unique(BAMM.mu.pgls.summary$Parameter)) {
hpd.Mu.top[[x]] = hdi(BAMM.mu.pgls.summary %>% filter(Parameter == x) %>% dplyr::select(Estimate))
}
hpd.Mu.top <- bind_rows(hpd.Mu.top) %>% `rownames<-`(c("Lower", "Upper")) %>% dplyr::select(-"(Intercept)")


saveRDS(hpd.Mu.top, 'data/hpd.Mu.top.rds')

hpd.DR.top %>% pander(split.table = Inf, digits = 3, caption = "DR HPD Interval")
DR HPD Interval
  SDi log(range.size.m2) bioclim4 residuals.PC1 PC1.LIG NPP
Lower -0.00193 -0.00844 -2.16e-05 -0.00323 -0.00389 -1.4e-06
Upper 0.00156 -0.00182 5.6e-05 0.00479 0.00523 1.53e-06
hpd.ND.top %>% pander(split.table = Inf, digits = 3, caption = "ND HPD Interval")
ND HPD Interval
  SDi log(range.size.m2) bioclim4 residuals.PC1 PC1.LIG NPP
Lower -7.68e-05 -0.00019 -9.72e-07 -9.09e-05 -0.00015 -4.27e-08
Upper 5.97e-05 9.11e-06 1.62e-06 0.000175 0.000117 4.28e-08
hpd.Lambda.top %>% pander(split.table = Inf, digits = 3, caption = "BAMM Speciation HPD Interval")
BAMM Speciation HPD Interval
  SDi log(range.size.m2) bioclim4 residuals.PC1 PC1.LIG NPP
Lower -0.00831 -0.00786 -0.000169 -0.0146 -0.0107 -5.56e-06
Upper 0.00474 0.0114 0.000124 0.0185 0.0084 2.42e-06
hpd.Mu.top %>% pander(split.table = Inf, digits = 3, caption = "BAMM Extinction HPD Interval")
BAMM Extinction HPD Interval
  SDi log(range.size.m2) bioclim4 residuals.PC1 PC1.LIG NPP
Lower -0.00788 -0.0165 -0.000159 -0.0216 -0.0161 -3.96e-06
Upper 0.0121 0.0237 0.000406 0.0367 0.0264 9.83e-06

From the HPD 95 % intervals and the intervals for the 95 % CI we can gain an estimate of how much variability there is within a model in comparison to variability between trees. To do this we calculate the width of the HPD interval relativel to the 95 % CI interval from the MCC.

Table S9: Model estimates from \(\lambda_{BAMM}\) and \(\mu_{BAMM}\) tip-rate estimates show high levels of variability between trees relative to estimates from \(\lambda_{DR}\) and \(\lambda_{ND}\). The interval ratio is calculated as the wisth of the HPD 95 % CI relative to the MCC 95 % CI. A value of 1 suggests that model estimates from 95 % of the trees fall within the MCC 95 % CI. Here, the large values for the estimates from BAMM tip-rates suggest that there is high variation in tip-rate estimates across trees that is not see within the DR and ND models.

hpd.DR.top2 <- hpd.DR.top %>% t() %>% as.data.frame() %>% rownames_to_column(var = "Parameter")
hpd.DR.top2$HPD.Interval <- hpd.DR.top2$Upper - hpd.DR.top2$Lower
DR.intervals <- left_join(hpd.DR.top2, MCC.DR.summary, by = "Parameter")
DR.intervals$MCC.Interval <- DR.intervals$UCI - DR.intervals$LCI
DR.intervals$DR.IntervalRatio <- hpd.DR.top2$HPD.Interval/DR.intervals$MCC.Interval

hpd.ND.top2 <- hpd.ND.top %>% t() %>% as.data.frame() %>% rownames_to_column(var = "Parameter")
hpd.ND.top2$HPD.Interval <- hpd.ND.top2$Upper - hpd.ND.top2$Lower
ND.intervals <- left_join(hpd.ND.top2, MCC.ND.summary, by = "Parameter")
ND.intervals$MCC.Interval <- ND.intervals$UCI - ND.intervals$LCI
ND.intervals$ND.IntervalRatio <- hpd.ND.top2$HPD.Interval/ND.intervals$MCC.Interval


hpd.Lambda.top2 <- hpd.Lambda.top %>% t() %>% as.data.frame() %>% rownames_to_column(var = "Parameter")
hpd.Lambda.top2$HPD.Interval <- hpd.Lambda.top2$Upper - hpd.Lambda.top2$Lower
Lambda.intervals <- left_join(hpd.Lambda.top2, MCC.lambda.summary, by = "Parameter")
Lambda.intervals$MCC.Interval <- Lambda.intervals$UCI - Lambda.intervals$LCI
Lambda.intervals$Lambda.IntervalRatio <- hpd.Lambda.top2$HPD.Interval/Lambda.intervals$MCC.Interval


hpd.Mu.top2 <- hpd.Mu.top %>% t() %>% as.data.frame() %>% rownames_to_column(var = "Parameter")
hpd.Mu.top2$HPD.Interval <- hpd.Mu.top2$Upper - hpd.Mu.top2$Lower
Mu.intervals <- left_join(hpd.Mu.top2, MCC.mu.summary, by = "Parameter")
Mu.intervals$MCC.Interval <- Mu.intervals$UCI - Mu.intervals$LCI
Mu.intervals$Mu.IntervalRatio <- hpd.Mu.top2$HPD.Interval/Mu.intervals$MCC.Interval


plyr::join_all(list(DR.intervals %>% select(Parameter, DR.IntervalRatio),
          ND.intervals %>% select(Parameter, ND.IntervalRatio),
          Lambda.intervals %>% select(Parameter, Lambda.IntervalRatio),
          Mu.intervals %>% select(Parameter, Mu.IntervalRatio)), by = "Parameter", type = "left") %>% `colnames<-`(c("Parameter", "λDR.Interval.Ratio", "λND.Interval.Ratio", "λBAMM.Interval.Ratio", "μBAMM.Interval.Ratio")) %>%
  pander(digits = 3, split.table = Inf) 
Parameter λDR.Interval.Ratio λND.Interval.Ratio λBAMM.Interval.Ratio μBAMM.Interval.Ratio
SDi 0.855 1.03 26.8 23.1
log(range.size.m2) 0.788 0.729 20.6 22.2
bioclim4 0.704 0.691 23.3 22.2
residuals.PC1 0.782 0.774 28.3 25.2
PC1.LIG 0.866 0.805 18.1 20.1
NPP 0.906 0.771 19.3 17.3

Subsetted analysis with Armenta data

Using the dataset from Armenta et al. (2008) we can conduct a subsetted analysis where dichromatism values are measured from spectrophotmetry and synthesised into a colour discriminability measure, which is a difference measure between the sexes. To make this dataset comprable with ours we use the absolute difference between the sexes. Therefore making the scale from monochromatism to dichromatism rather than female colouration to male colouration.

#Read in Armenta data
Armenta.data <- read.csv('data/Armenta_2008.csv', stringsAsFactors = F)
MCC.BAMM.df <- readRDS('data/MCC.BAMM.df.rds')
#A couple of rows have "-", remove from dataset
Armenta.data <- Armenta.data %>% dplyr::select(binomial, Colour.discriminability) %>% mutate(Colour.discriminability = replace(Colour.discriminability, Colour.discriminability == "-", "NA")) %>% filter(Colour.discriminability != "NA")
Armenta.data$Colour.discriminability <- Armenta.data$Colour.discriminability %>% as.numeric()

MCC.Armenta.model.data <- inner_join(restricted.data %>% dplyr::select(TipLabel, binomial, SDi, range.size.m2, bioclim4, residuals.PC1, PC1.LIG, NPP, MCC.DR, MCC.ND), Armenta.data, by = "binomial")

MCC.Armenta.model.data <- inner_join(MCC.Armenta.model.data, MCC.BAMM.df, by = "TipLabel") %>% filter(Colour.discriminability != "NA")

#Plot the correlation
p1 <- MCC.Armenta.model.data %>% ggplot(aes(x = SDi, y = Colour.discriminability))+
  geom_point()+
  geom_smooth(method = 'loess')+
  theme_minimal()

grid.newpage()
ggExtra::ggMarginal(
  p = p1,
  type = 'density',
  margins = 'both',
  size = 5,
  colour = 'black',
  fill = 'gray'
)

Figure S9: There is a correlation between the RGB measures and the Spectophotometry measures, The RGB measures seem to be more noisy around the lower values.

R-squared and correlation for relationship between RGB and Spectrophotometry Data:

data_frame(
R2 = summary(lm(SDi ~ Colour.discriminability,
   data = MCC.Armenta.model.data))$r.squared,
r = cor(MCC.Armenta.model.data$SDi, MCC.Armenta.model.data$Colour.discriminability)) %>% pander()
R2 r
0.6166 0.7853

Run PGLS model for the data

#Prune tree
pruned.MCC.Armenta.tree <- drop.tip(MCC.passerine,MCC.passerine$tip.label[-match(MCC.Armenta.model.data$TipLabel, MCC.passerine$tip.label)])

#Set rownames to match tree
rownames(MCC.Armenta.model.data) <- MCC.Armenta.model.data$TipLabel

#Run a corPagel model to estimate lambda for DR
MCC.DR.Armenta.global <- gls(MCC.DR ~ Colour.discriminability
                         + log(range.size.m2)
                         + bioclim4 #Seasonal variation
                         + residuals.PC1 #Spatial variation
                         + PC1.LIG #Long-term climate variation
                         + NPP
                         + Colour.discriminability*log(range.size.m2)
                         + Colour.discriminability*bioclim4
                         + Colour.discriminability*residuals.PC1
                         + Colour.discriminability*PC1.LIG
                         + Colour.discriminability*NPP,
                correlation = corPagel(1, phy = pruned.MCC.Armenta.tree, fixed = FALSE), 
                data = MCC.Armenta.model.data, 
                method = "REML")

#Run a corPagel model to estimate lambda for DR
MCC.ND.Armenta.global <- gls(MCC.ND ~ Colour.discriminability
                         + log(range.size.m2)
                         + bioclim4 #Seasonal variation
                         + residuals.PC1 #Spatial variation
                         + PC1.LIG #Long-term climate variation
                         + NPP
                         + Colour.discriminability*log(range.size.m2)
                         + Colour.discriminability*bioclim4
                         + Colour.discriminability*residuals.PC1
                         + Colour.discriminability*PC1.LIG
                         + Colour.discriminability*NPP,
                correlation = corPagel(1, phy = pruned.MCC.Armenta.tree, fixed = TRUE), #lambda = 1
                data = MCC.Armenta.model.data, 
                method = "REML")

#Run a corPagel model to estimate lambda for DR
MCC.Lambda.Armenta.global <- gls(mean.lambda ~ Colour.discriminability
                         + log(range.size.m2)
                         + bioclim4 #Seasonal variation
                         + residuals.PC1 #Spatial variation
                         + PC1.LIG #Long-term climate variation
                         + NPP
                         + Colour.discriminability*log(range.size.m2)
                         + Colour.discriminability*bioclim4
                         + Colour.discriminability*residuals.PC1
                         + Colour.discriminability*PC1.LIG
                         + Colour.discriminability*NPP,
                weights = ~ sqrt(var.lambda),
                correlation = corPagel(1, phy = pruned.MCC.Armenta.tree, fixed = TRUE), 
                data = MCC.Armenta.model.data, 
                method = "REML")

#Run a corPagel model to estimate lambda for DR
MCC.Extinction.Armenta.global <- gls(mean.mu ~ Colour.discriminability
                         + log(range.size.m2)
                         + bioclim4 #Seasonal variation
                         + residuals.PC1 #Spatial variation
                         + PC1.LIG #Long-term climate variation
                         + NPP
                         + Colour.discriminability*log(range.size.m2)
                         + Colour.discriminability*bioclim4
                         + Colour.discriminability*residuals.PC1
                         + Colour.discriminability*PC1.LIG
                         + Colour.discriminability*NPP,
                weights = ~ sqrt(var.mu),
                correlation = corPagel(1, phy = pruned.MCC.Armenta.tree, fixed = TRUE), 
                data = MCC.Armenta.model.data, 
                method = "REML")
#Set up cluster
cores<-8
clusterType <- if(length(find.package("snow", quiet = TRUE))) "SOCK" else "PSOCK"
clust <- try(makeCluster(getOption("cl.cores", cores), type = clusterType))
myclust<-clust

#Export data and packages to cluster
clusterExport(myclust, c("MCC.Armenta.model.data"), envir=environment())
clusterExport(myclust, c("pruned.MCC.Armenta.tree"), envir=environment())
clusterEvalQ(myclust, library(nlme))
clusterEvalQ(myclust, library(ape))
clusterEvalQ(myclust, library(MuMIn))

#Dredged models:

Armenta.dredged.ND.model <- pdredge(MCC.ND.Armenta.global, fixed = c("Colour.discriminability", "log(range.size.m2)", "bioclim4", "residuals.PC1", "PC1.LIG", "NPP"), cluster=myclust)

saveRDS(Armenta.dredged.ND.model, 'data/Armenta.dredged.ND.model.rds')

Armenta.dredged.DR.model <- pdredge(MCC.DR.Armenta.global, fixed = c("Colour.discriminability", "log(range.size.m2)", "bioclim4", "residuals.PC1", "PC1.LIG", "NPP"), cluster=myclust)

saveRDS(Armenta.dredged.DR.model, 'data/Armenta.dredged.DR.model.rds')

Armenta.dredged.spec.model <- pdredge(MCC.Lambda.Armenta.global, fixed = c("Colour.discriminability", "log(range.size.m2)", "bioclim4", "residuals.PC1", "PC1.LIG", "NPP"), cluster=myclust)

saveRDS(Armenta.dredged.spec.model, 'data/Armenta.dredged.spec.model.rds')

Armenta.dredged.extinct.model <- pdredge(MCC.Extinction.Armenta.global, fixed = c("Colour.discriminability", "log(range.size.m2)", "bioclim4", "residuals.PC1", "PC1.LIG", "NPP"), cluster=myclust)

saveRDS(Armenta.dredged.extinct.model, 'data/Armenta.dredged.extinct.model.rds')

Table S10: The dredged models for Models using a subset dataset where sexual dichromatism was measured using spectrophotometry, all four show the top model is one with no interactions , with \(\delta AICc > 4\).

Armenta.dredged.ND.model <- readRDS("data/Armenta.dredged.ND.model.rds")
Armenta.dredged.DR.model <- readRDS("data/Armenta.dredged.DR.model.rds")
Armenta.dredged.spec.model <- readRDS("data/Armenta.dredged.spec.model.rds")
Armenta.dredged.extinct.model <- readRDS("data/Armenta.dredged.extinct.model.rds")


kable(Armenta.dredged.ND.model, "html", caption = "ND Spectrophotometry Dichromatism Dredge Table") %>%
  kable_styling() %>%
  scroll_box(width = "100%", height = "500px")
ND Spectrophotometry Dichromatism Dredge Table
(Intercept) bioclim4 Colour.discriminability log(range.size.m2) NPP PC1.LIG residuals.PC1 bioclim4:Colour.discriminability Colour.discriminability:log(range.size.m2) Colour.discriminability:NPP Colour.discriminability:PC1.LIG Colour.discriminability:residuals.PC1 df logLik AICc delta weight
0 0.2077881 6.9e-06 0.0000828 -0.0000378 5e-07 -0.0008841 0.0010532 NA NA NA NA NA 8 1033.2946 -2050.337 0.00000 9.972312e-01
8 0.2128993 6.2e-06 0.0003666 -0.0002044 4e-07 -0.0008428 0.0010460 NA NA NA -0.0009029 NA 9 1027.8860 -2037.457 12.88059 1.591519e-03
2 0.2060699 7.3e-06 -0.0081861 0.0000182 5e-07 -0.0008841 0.0010718 NA 0.0002887 NA NA NA 9 1027.0481 -2035.781 14.55654 6.884668e-04
16 0.2083870 6.9e-06 0.0001261 -0.0000596 5e-07 -0.0008744 0.0010533 NA NA NA NA -0.0001539 9 1026.6890 -2035.063 15.27474 4.807590e-04
1 0.2072627 7.1e-06 -0.0002087 -0.0000199 5e-07 -0.0008882 0.0010535 6.00e-07 NA NA NA NA 9 1021.9461 -2025.577 24.76044 4.189235e-06
10 0.2094193 7.0e-06 -0.0259111 -0.0000910 5e-07 -0.0008267 0.0011025 NA 0.0009215 NA -0.0012568 NA 10 1022.3412 -2024.296 26.04102 2.208311e-06
24 0.2127267 6.2e-06 0.0003401 -0.0001949 4e-07 -0.0008524 0.0010447 NA NA NA -0.0010503 0.0002588 10 1021.4594 -2022.533 27.80461 9.143286e-07
4 0.2077525 7.1e-06 0.0008792 -0.0000384 5e-07 -0.0009151 0.0010575 NA NA -1e-07 NA NA 9 1019.5502 -2020.785 29.55223 3.816023e-07
18 0.2066684 7.2e-06 -0.0081477 -0.0000036 5e-07 -0.0008744 0.0010719 NA 0.0002889 NA NA -0.0001540 10 1020.4432 -2020.500 29.83700 3.309587e-07
9 0.2095915 7.9e-06 -0.0054207 -0.0000803 4e-07 -0.0008660 0.0010408 1.38e-05 NA NA -0.0023002 NA 10 1018.8039 -2017.222 33.11555 6.424583e-08
3 0.2062134 7.2e-06 -0.0088719 0.0000131 5e-07 -0.0008816 0.0010736 -4.00e-07 0.0003190 NA NA NA 10 1015.8364 -2011.287 39.05055 3.304277e-09
17 0.2074141 7.1e-06 -0.0005561 -0.0000272 5e-07 -0.0008795 0.0010540 1.60e-06 NA NA NA -0.0002296 10 1015.4839 -2010.582 39.75568 2.322516e-09
26 0.2084783 7.2e-06 -0.0307006 -0.0000533 5e-07 -0.0008411 0.0011104 NA 0.0010878 NA -0.0015881 0.0004694 11 1016.1516 -2009.839 40.49816 1.602260e-09
12 0.2133899 6.8e-06 0.0028641 -0.0002261 4e-07 -0.0009338 0.0010583 NA NA -3e-07 -0.0010090 NA 10 1014.3466 -2008.307 42.03013 7.448502e-10
6 0.2060643 7.3e-06 -0.0082262 0.0000185 5e-07 -0.0008837 0.0010718 NA 0.0002898 0e+00 NA NA 10 1013.4370 -2006.488 43.84935 2.999366e-10
20 0.2083692 7.1e-06 0.0009782 -0.0000611 5e-07 -0.0009072 0.0010578 NA NA -1e-07 NA -0.0001591 10 1012.9494 -2005.513 44.82470 1.841768e-10
11 0.2078948 8.2e-06 -0.0201410 -0.0000266 4e-07 -0.0008543 0.0010742 1.24e-05 0.0005364 NA -0.0023673 NA 11 1012.8327 -2003.202 47.13593 5.799078e-11
25 0.2095645 7.9e-06 -0.0053564 -0.0000780 4e-07 -0.0008695 0.0010404 1.36e-05 NA NA -0.0023418 0.0001047 11 1012.3028 -2002.142 48.19581 3.413571e-11
19 0.2065718 7.2e-06 -0.0071738 -0.0000002 5e-07 -0.0008763 0.0010694 6.00e-07 0.0002463 NA NA -0.0001812 11 1009.3819 -1996.300 54.03762 1.839399e-12
5 0.2074923 7.2e-06 0.0006301 -0.0000293 5e-07 -0.0009132 0.0010571 3.00e-07 NA -1e-07 NA NA 10 1008.2734 -1996.161 54.17659 1.715926e-12
14 0.2096468 7.1e-06 -0.0243990 -0.0000994 4e-07 -0.0008428 0.0011022 NA 0.0008833 0e+00 -0.0012602 NA 11 1008.7339 -1995.004 55.33369 9.621384e-13
28 0.2132308 6.8e-06 0.0030279 -0.0002169 4e-07 -0.0009518 0.0010578 NA NA -3e-07 -0.0011857 0.0002957 11 1007.9537 -1993.443 56.89410 4.409596e-13
22 0.2067012 7.2e-06 -0.0079237 -0.0000048 5e-07 -0.0008768 0.0010718 NA 0.0002832 0e+00 NA -0.0001544 11 1006.8344 -1991.205 59.13254 1.439887e-13
27 0.2074778 8.2e-06 -0.0230409 -0.0000101 4e-07 -0.0008604 0.0010801 1.17e-05 0.0006477 NA -0.0024811 0.0002517 12 1006.4510 -1988.353 61.98466 3.459381e-14
13 0.2096921 7.9e-06 -0.0050367 -0.0000845 4e-07 -0.0008772 0.0010424 1.36e-05 NA 0e+00 -0.0022975 NA 11 1005.1167 -1987.770 62.56792 2.584311e-14
7 0.2062406 7.2e-06 -0.0087170 0.0000121 5e-07 -0.0008834 0.0010736 -4.00e-07 0.0003153 0e+00 NA NA 11 1002.2404 -1982.017 68.32055 1.456052e-15
21 0.2075602 7.2e-06 0.0000142 -0.0000332 5e-07 -0.0008964 0.0010563 1.30e-06 NA -1e-07 NA -0.0002203 11 1001.8148 -1981.166 69.17175 9.513509e-16
30 0.2086525 7.3e-06 -0.0295459 -0.0000597 5e-07 -0.0008532 0.0011102 NA 0.0010585 0e+00 -0.0015897 0.0004681 12 1002.5432 -1980.537 69.80031 6.947876e-16
15 0.2073335 8.1e-06 -0.0234976 -0.0000057 4e-07 -0.0008176 0.0010741 1.27e-05 0.0006154 1e-07 -0.0023857 NA 12 999.2588 -1973.968 76.36914 2.602837e-17
29 0.2096993 7.9e-06 -0.0048284 -0.0000836 4e-07 -0.0008852 0.0010425 1.33e-05 NA 0e+00 -0.0023414 0.0001133 12 998.6289 -1972.708 77.62895 1.386383e-17
23 0.2065684 7.2e-06 -0.0071929 -0.0000001 5e-07 -0.0008760 0.0010694 6.00e-07 0.0002468 0e+00 NA -0.0001813 12 995.7873 -1967.025 83.31204 8.087520e-19
31 0.2069328 8.1e-06 -0.0262977 0.0000102 4e-07 -0.0008246 0.0010799 1.20e-05 0.0007241 1e-07 -0.0024983 0.0002501 13 992.8768 -1959.112 91.22581 1.546546e-20
kable(Armenta.dredged.DR.model, "html", caption = "DR Spectrophotometry Dichromatism Dredge Table") %>%
  kable_styling() %>%
  scroll_box(width = "100%", height = "500px")
DR Spectrophotometry Dichromatism Dredge Table
(Intercept) bioclim4 Colour.discriminability log(range.size.m2) NPP PC1.LIG residuals.PC1 bioclim4:Colour.discriminability Colour.discriminability:log(range.size.m2) Colour.discriminability:NPP Colour.discriminability:PC1.LIG Colour.discriminability:residuals.PC1 df logLik AICc delta weight
0 -2.677136 0.0001789 0.0070708 0.0090337 7.8e-06 -0.0315655 0.0156001 NA NA NA NA NA 9 -512.3837 1043.083 0.000000 9.479871e-01
8 -2.630812 0.0001604 0.0092806 0.0078425 7.2e-06 -0.0288585 0.0166944 NA NA NA -0.0169326 NA 10 -514.7327 1049.851 6.768670 3.213693e-02
2 -2.669309 0.0001769 0.0667921 0.0087960 7.7e-06 -0.0314417 0.0155465 NA -0.0020932 NA NA NA 10 -515.8406 1052.067 8.984518 1.061303e-02
16 -2.675941 0.0001777 0.0078042 0.0089905 7.8e-06 -0.0313098 0.0157537 NA NA NA NA -0.0017242 10 -516.0987 1052.583 9.500619 8.199156e-03
10 -2.653060 0.0001640 -0.2224726 0.0085957 7.3e-06 -0.0288185 0.0171615 NA 0.0081372 NA -0.0199731 NA 11 -517.9245 1058.313 15.230367 4.672726e-04
24 -2.624359 0.0001602 0.0064975 0.0077898 6.9e-06 -0.0292878 0.0163567 NA NA NA -0.0213744 0.0079966 11 -518.0589 1058.582 15.499006 4.085413e-04
18 -2.668225 0.0001757 0.0668268 0.0087564 7.8e-06 -0.0311887 0.0157002 NA -0.0020689 NA NA -0.0017135 11 -519.5546 1061.573 18.490405 9.155078e-05
1 -2.671939 0.0001754 0.0136694 0.0088999 7.8e-06 -0.0312585 0.0157116 -0.0000153 NA NA NA NA 10 -520.8980 1062.182 19.099199 6.752492e-05
9 -2.642496 0.0001826 -0.0714208 0.0081999 7.2e-06 -0.0295284 0.0166390 0.0001929 NA NA -0.0366791 NA 11 -521.5732 1065.610 22.527776 1.216068e-05
4 -2.673868 0.0001865 0.0487209 0.0088595 7.7e-06 -0.0331665 0.0160592 NA NA -4.6e-06 NA NA 10 -523.1346 1066.655 23.572501 7.212717e-06
26 -2.654133 0.0001654 -0.3246079 0.0088415 7.0e-06 -0.0293710 0.0169056 NA 0.0115950 NA -0.0270764 0.0104643 12 -521.0596 1066.669 23.585830 7.164808e-06
3 -2.669100 0.0001751 0.0435357 0.0088140 7.7e-06 -0.0312622 0.0156661 -0.0000116 -0.0011029 NA NA NA 11 -524.1891 1070.842 27.759554 8.889798e-07
17 -2.672419 0.0001753 0.0125561 0.0089088 7.8e-06 -0.0311560 0.0157961 -0.0000115 NA NA NA -0.0011657 11 -524.5318 1071.528 28.444956 6.310427e-07
12 -2.618510 0.0001699 0.0766749 0.0073776 6.8e-06 -0.0310209 0.0175985 NA NA -7.4e-06 -0.0195276 NA 11 -525.1224 1072.709 29.626063 3.496106e-07
11 -2.643111 0.0001824 -0.0776207 0.0082299 7.2e-06 -0.0295154 0.0166688 0.0001919 0.0002323 NA -0.0366674 NA 12 -524.8734 1074.296 31.213509 1.580793e-07
25 -2.636482 0.0001821 -0.0727080 0.0081501 6.9e-06 -0.0299139 0.0163324 0.0001898 NA NA -0.0404677 0.0073866 12 -524.9401 1074.429 31.346776 1.478892e-07
6 -2.638321 0.0001829 0.3290294 0.0077113 7.3e-06 -0.0335589 0.0160458 NA -0.0090012 -7.2e-06 NA NA 11 -526.2544 1074.973 31.890077 1.127094e-07
20 -2.672935 0.0001854 0.0488811 0.0088282 7.7e-06 -0.0329441 0.0161800 NA NA -4.6e-06 NA -0.0013734 11 -526.8525 1076.169 33.086302 6.197310e-08
19 -2.668293 0.0001750 0.0571451 0.0087749 7.7e-06 -0.0311417 0.0157386 -0.0000048 -0.0016595 NA NA -0.0014888 12 -527.7809 1080.111 37.028433 8.633371e-09
28 -2.607544 0.0001714 0.0853530 0.0072208 6.4e-06 -0.0320005 0.0173071 NA NA -8.8e-06 -0.0259988 0.0107877 12 -528.2316 1081.012 37.929758 5.501233e-09
14 -2.622604 0.0001699 0.0367068 0.0075259 6.9e-06 -0.0309141 0.0176481 NA 0.0013029 -7.1e-06 -0.0199055 NA 12 -528.3498 1081.249 38.166275 4.887662e-09
27 -2.644834 0.0001822 -0.1679306 0.0084555 7.0e-06 -0.0298937 0.0165119 0.0001781 0.0035063 NA -0.0410063 0.0081624 13 -528.1672 1082.976 39.893793 2.060513e-09
22 -2.637873 0.0001820 0.3270864 0.0076981 7.3e-06 -0.0333773 0.0161453 NA -0.0089347 -7.2e-06 NA -0.0010874 12 -529.9733 1084.496 41.413285 9.638776e-10
5 -2.656263 0.0001785 0.0857028 0.0083648 7.4e-06 -0.0328347 0.0165795 -0.0000478 NA -6.4e-06 NA NA 11 -531.4166 1085.297 42.214478 6.457213e-10
13 -2.635478 0.0001841 -0.0303738 0.0079455 7.0e-06 -0.0304341 0.0170601 0.0001683 NA -3.4e-06 -0.0353444 NA 12 -532.3431 1089.235 46.152785 9.012648e-11
30 -2.620812 0.0001720 -0.0499593 0.0076922 6.5e-06 -0.0317146 0.0174172 NA 0.0044218 -7.8e-06 -0.0276394 0.0114010 13 -531.4062 1089.454 46.371621 8.078539e-11
7 -2.635314 0.0001785 0.2873035 0.0076666 7.2e-06 -0.0332479 0.0163926 -0.0000308 -0.0068943 -7.8e-06 NA NA 12 -534.5375 1093.624 50.541586 1.004237e-11
21 -2.654840 0.0001782 0.0907471 0.0083515 7.3e-06 -0.0330081 0.0165591 -0.0000543 NA -6.7e-06 NA 0.0014471 12 -535.0039 1094.557 51.474387 6.299140e-12
15 -2.627129 0.0001839 0.0527643 0.0076688 6.9e-06 -0.0306250 0.0169889 0.0001727 -0.0027968 -4.0e-06 -0.0349658 NA 13 -535.5408 1097.724 54.640887 1.293256e-12
29 -2.624974 0.0001842 -0.0138454 0.0077690 6.6e-06 -0.0313081 0.0168654 0.0001537 NA -4.9e-06 -0.0393974 0.0090385 13 -535.5861 1097.814 54.731532 1.235951e-12
23 -2.635166 0.0001782 0.2848941 0.0076832 7.2e-06 -0.0333049 0.0164014 -0.0000341 -0.0067339 -7.9e-06 NA 0.0006111 13 -538.1141 1102.870 59.787439 9.865632e-14
31 -2.625387 0.0001840 -0.0182030 0.0077919 6.6e-06 -0.0312931 0.0168860 0.0001531 0.0001569 -4.9e-06 -0.0393997 0.0090593 14 -538.7716 1106.285 63.202639 1.788643e-14
kable(Armenta.dredged.spec.model, "html", caption = "BAMM Speciation Spectrophotometry Dichromatism Dredge Table") %>%
  kable_styling() %>%
  scroll_box(width = "100%", height = "500px")
BAMM Speciation Spectrophotometry Dichromatism Dredge Table
(Intercept) bioclim4 Colour.discriminability log(range.size.m2) NPP PC1.LIG residuals.PC1 bioclim4:Colour.discriminability Colour.discriminability:log(range.size.m2) Colour.discriminability:NPP Colour.discriminability:PC1.LIG Colour.discriminability:residuals.PC1 df logLik AICc delta weight
0 -2.482091 -2.59e-05 -0.0005285 0.0018448 -1.6e-06 -0.0042436 0.0024890 NA NA NA NA NA 8 181.4767 -346.7017 0.00000 9.919991e-01
2 -2.488239 -2.42e-05 -0.0427556 0.0020142 -1.6e-06 -0.0042772 0.0025602 NA 0.0014676 NA NA NA 9 176.7474 -335.1796 11.52204 3.122707e-03
8 -2.477350 -2.71e-05 -0.0001200 0.0017133 -1.7e-06 -0.0041157 0.0025958 NA NA NA -0.0011826 NA 9 176.5780 -334.8407 11.86100 2.635893e-03
16 -2.478790 -2.64e-05 -0.0000802 0.0017377 -1.6e-06 -0.0040974 0.0025735 NA NA NA NA -0.0009394 9 176.3953 -334.4753 12.22638 2.195774e-03
1 -2.475060 -2.79e-05 0.0040917 0.0016355 -1.7e-06 -0.0041107 0.0025852 -1.04e-05 NA NA NA NA 9 171.7160 -325.1168 21.58489 2.038976e-05
10 -2.484288 -2.53e-05 -0.0750059 0.0019038 -1.6e-06 -0.0040674 0.0028130 NA 0.0026147 NA -0.0021843 NA 10 172.1264 -323.8668 22.83488 1.091389e-05
18 -2.484933 -2.47e-05 -0.0422043 0.0019071 -1.6e-06 -0.0041314 0.0026442 NA 0.0014640 NA NA -0.0009365 10 171.6660 -322.9461 23.75558 6.887368e-06
24 -2.476734 -2.71e-05 0.0000344 0.0016866 -1.7e-06 -0.0040638 0.0026159 NA NA NA -0.0008273 -0.0005807 10 171.5901 -322.7942 23.90753 6.383472e-06
4 -2.482182 -2.61e-05 -0.0016774 0.0018484 -1.6e-06 -0.0042020 0.0024819 NA NA 1e-07 NA NA 9 169.2317 -320.1481 26.55354 1.700135e-06
3 -2.481241 -2.64e-05 -0.0781279 0.0017918 -1.6e-06 -0.0040585 0.0028184 -1.98e-05 0.0030037 NA NA NA 10 167.3597 -314.3334 32.36824 9.286167e-08
9 -2.475019 -2.79e-05 0.0034649 0.0016366 -1.7e-06 -0.0041007 0.0025958 -8.70e-06 NA NA -0.0002856 NA 10 167.1645 -313.9430 32.75864 7.639465e-08
17 -2.474660 -2.78e-05 0.0034158 0.0016188 -1.7e-06 -0.0040574 0.0026121 -8.30e-06 NA NA NA -0.0005125 10 166.6862 -312.9865 33.71517 4.735391e-08
26 -2.484329 -2.53e-05 -0.0752360 0.0019053 -1.6e-06 -0.0040689 0.0028130 NA 0.0026225 NA -0.0021992 0.0000193 11 167.1639 -311.8639 34.83783 2.701296e-08
6 -2.491163 -2.51e-05 -0.0658585 0.0021010 -1.6e-06 -0.0040506 0.0025469 NA 0.0020406 7e-07 NA NA 10 164.6871 -308.9882 37.71345 6.414165e-09
12 -2.477129 -2.70e-05 0.0008182 0.0017063 -1.7e-06 -0.0041453 0.0026049 NA NA -1e-07 -0.0012197 NA 10 164.3681 -308.3503 38.35143 4.662337e-09
20 -2.478881 -2.67e-05 -0.0012593 0.0017414 -1.6e-06 -0.0040546 0.0025662 NA NA 1e-07 NA -0.0009401 10 164.1511 -307.9162 38.78547 3.752787e-09
11 -2.481302 -2.63e-05 -0.0819535 0.0017988 -1.6e-06 -0.0040313 0.0028522 -1.59e-05 0.0030845 NA -0.0007358 NA 11 162.8293 -303.1947 43.50702 3.540640e-10
19 -2.481434 -2.64e-05 -0.0793184 0.0017982 -1.6e-06 -0.0040698 0.0028161 -2.04e-05 0.0030529 NA NA 0.0001175 11 162.3563 -302.2485 44.45314 2.206145e-10
25 -2.474659 -2.78e-05 0.0033284 0.0016192 -1.7e-06 -0.0040567 0.0026133 -8.10e-06 NA NA -0.0000445 -0.0005046 11 162.1775 -301.8910 44.81064 1.845031e-10
5 -2.474008 -2.75e-05 0.0074117 0.0016021 -1.7e-06 -0.0041954 0.0026138 -1.16e-05 NA -3e-07 NA NA 10 159.5612 -298.7363 47.96534 3.810396e-11
14 -2.486942 -2.61e-05 -0.0945190 0.0019825 -1.6e-06 -0.0038735 0.0027955 NA 0.0030906 6e-07 -0.0021341 NA 11 160.0554 -297.6469 49.05480 2.210023e-11
22 -2.487860 -2.57e-05 -0.0653996 0.0019940 -1.6e-06 -0.0039035 0.0026310 NA 0.0020393 7e-07 NA -0.0009389 11 159.6068 -296.7496 49.95208 1.411094e-11
28 -2.476677 -2.70e-05 0.0003005 0.0016849 -1.7e-06 -0.0040728 0.0026183 NA NA 0e+00 -0.0008419 -0.0005743 11 159.3941 -296.3242 50.37745 1.140740e-11
27 -2.481836 -2.63e-05 -0.0860582 0.0018177 -1.6e-06 -0.0040555 0.0028538 -1.66e-05 0.0032358 NA -0.0009089 0.0003162 12 157.8956 -291.2420 55.45971 8.986463e-13
7 -2.483006 -2.68e-05 -0.0885776 0.0018448 -1.6e-06 -0.0039510 0.0028016 -1.91e-05 0.0032375 4e-07 NA NA 11 155.2786 -288.0933 58.60841 1.861469e-13
13 -2.473997 -2.75e-05 0.0067813 0.0016039 -1.7e-06 -0.0041846 0.0026225 -1.01e-05 NA -3e-07 -0.0002506 NA 11 155.0117 -287.5594 59.14227 1.425376e-13
21 -2.473936 -2.75e-05 0.0059245 0.0015963 -1.7e-06 -0.0041256 0.0026299 -9.40e-06 NA -2e-07 NA -0.0004528 11 154.5518 -286.6396 60.06208 8.999034e-14
30 -2.486908 -2.61e-05 -0.0943380 0.0019813 -1.6e-06 -0.0038720 0.0027955 NA 0.0030843 6e-07 -0.0021212 -0.0000167 12 155.0948 -285.6402 61.06147 5.459869e-14
15 -2.483310 -2.67e-05 -0.0942672 0.0018598 -1.6e-06 -0.0039062 0.0028373 -1.45e-05 0.0033594 4e-07 -0.0008255 NA 12 150.7592 -276.9690 69.73265 7.149220e-16
23 -2.483076 -2.68e-05 -0.0890223 0.0018472 -1.6e-06 -0.0039590 0.0028007 -1.94e-05 0.0032594 4e-07 NA 0.0000624 12 150.2835 -276.0176 70.68404 4.442896e-16
29 -2.473936 -2.75e-05 0.0058337 0.0015967 -1.7e-06 -0.0041249 0.0026312 -9.20e-06 NA -2e-07 -0.0000470 -0.0004445 12 150.0439 -275.5386 71.16310 3.496553e-16
31 -2.483668 -2.67e-05 -0.0971660 0.0018729 -1.6e-06 -0.0039327 0.0028394 -1.52e-05 0.0034748 4e-07 -0.0009683 0.0002686 13 145.8290 -265.0161 81.68558 1.814318e-18
kable(Armenta.dredged.extinct.model, "html", caption = "BAMM Extinction Spectrophotometry Dichromatism Dredge Table") %>%
  kable_styling() %>%
  scroll_box(width = "100%", height = "500px")
BAMM Extinction Spectrophotometry Dichromatism Dredge Table
(Intercept) bioclim4 Colour.discriminability log(range.size.m2) NPP PC1.LIG residuals.PC1 bioclim4:Colour.discriminability Colour.discriminability:log(range.size.m2) Colour.discriminability:NPP Colour.discriminability:PC1.LIG Colour.discriminability:residuals.PC1 df logLik AICc delta weight
0 -5.638235 -2.60e-06 -0.0007063 0.0098190 -7.0e-07 -0.0158136 0.0047533 NA NA NA NA NA 8 -295.8250 607.9017 0.000000 9.760815e-01
8 -5.603206 -9.50e-06 0.0013768 0.0086377 -1.0e-06 -0.0153132 0.0048284 NA NA NA -0.0069542 NA 9 -299.3662 617.0476 9.145821 1.008083e-02
2 -5.660518 3.00e-06 -0.1231044 0.0104727 -5.0e-07 -0.0158420 0.0050480 NA 0.0042674 NA NA NA 9 -299.6673 617.6499 9.748116 7.459498e-03
16 -5.622155 -4.50e-06 0.0006290 0.0091959 -8.0e-07 -0.0155208 0.0048853 NA NA NA NA -0.0039732 9 -299.8737 618.0626 10.160901 6.068405e-03
10 -5.634016 -7.00e-07 -0.2647213 0.0094634 -8.0e-07 -0.0151204 0.0055100 NA 0.0093146 NA -0.0104934 NA 10 -302.5940 625.5740 17.672233 1.419088e-04
1 -5.607926 -1.09e-05 0.0181308 0.0087106 -8.0e-07 -0.0155487 0.0048392 -4.15e-05 NA NA NA NA 9 -304.4774 627.2700 19.368277 6.077406e-05
24 -5.601118 -9.30e-06 0.0016660 0.0085321 -1.0e-06 -0.0152558 0.0048753 NA NA NA -0.0059789 -0.0017298 10 -303.5448 627.4756 19.573860 5.483734e-05
18 -5.643865 9.00e-07 -0.1162808 0.0098369 -6.0e-07 -0.0155557 0.0051632 NA 0.0040748 NA NA -0.0038665 10 -303.7310 627.8479 19.946169 4.552297e-05
4 -5.638385 1.60e-06 0.0142457 0.0097934 -7.0e-07 -0.0164385 0.0048215 NA NA -1.6e-06 NA NA 9 -307.2012 632.7177 24.815946 3.988153e-06
3 -5.634989 -4.40e-06 -0.2351668 0.0093492 -6.0e-07 -0.0154226 0.0055422 -7.09e-05 0.0092968 NA NA NA 10 -307.7169 635.8199 27.918123 8.455567e-07
26 -5.634738 -6.00e-07 -0.2680032 0.0094910 -8.0e-07 -0.0151276 0.0055105 NA 0.0094278 NA -0.0106989 0.0002882 11 -306.7717 636.0073 28.105536 7.699216e-07
9 -5.601830 -1.01e-05 0.0038916 0.0085790 -1.0e-06 -0.0153175 0.0048344 -5.90e-06 NA NA -0.0063678 NA 10 -307.9890 636.3640 28.462222 6.441591e-07
17 -5.605057 -1.01e-05 0.0144751 0.0085837 -9.0e-07 -0.0154279 0.0049016 -3.16e-05 NA NA NA -0.0024957 10 -308.6217 637.6294 29.727669 3.421408e-07
12 -5.597342 -2.80e-06 0.0298211 0.0083824 -1.2e-06 -0.0163988 0.0049698 NA NA -3.1e-06 -0.0081742 NA 10 -310.5363 641.4586 33.556851 5.043227e-08
6 -5.658376 3.70e-06 -0.1065804 0.0104008 -5.0e-07 -0.0160276 0.0050397 NA 0.0038485 -5.0e-07 NA NA 10 -310.9671 642.3201 34.418384 3.278151e-08
20 -5.621991 0.00e+00 0.0167326 0.0091558 -8.0e-07 -0.0161867 0.0049613 NA NA -1.8e-06 NA -0.0040535 10 -311.2374 642.8608 34.959076 2.501608e-08
11 -5.629948 -2.80e-06 -0.2763034 0.0092454 -7.0e-07 -0.0151246 0.0056014 -2.96e-05 0.0101592 NA -0.0078889 NA 11 -311.1065 644.6769 36.775141 1.008941e-08
19 -5.633651 -4.40e-06 -0.2294802 0.0093043 -6.0e-07 -0.0153988 0.0055382 -6.79e-05 0.0090580 NA NA -0.0005592 11 -311.9033 646.2705 38.368752 4.547968e-09
25 -5.600503 -9.60e-06 0.0028601 0.0085062 -1.0e-06 -0.0152591 0.0048772 -2.80e-06 NA NA -0.0057187 -0.0016952 11 -312.1582 646.7804 38.878672 3.524436e-09
14 -5.629632 7.00e-07 -0.2337836 0.0093153 -9.0e-07 -0.0154751 0.0054976 NA 0.0085422 -1.0e-06 -0.0105806 NA 11 -313.8788 650.2216 42.319887 6.307240e-10
5 -5.597683 -4.80e-06 0.0563682 0.0082702 -1.0e-06 -0.0167800 0.0050137 -5.60e-05 NA -3.5e-06 NA NA 10 -315.5658 651.5177 43.615917 3.299210e-10
28 -5.595913 -2.90e-06 0.0291288 0.0083088 -1.2e-06 -0.0163194 0.0050015 NA NA -3.0e-06 -0.0073794 -0.0013401 11 -314.7274 651.9188 44.017031 2.699661e-10
22 -5.640481 1.90e-06 -0.0917470 0.0097224 -6.0e-07 -0.0158266 0.0051523 NA 0.0034525 -7.0e-07 NA -0.0039161 11 -315.0223 652.5086 44.606904 2.010110e-10
27 -5.631914 -2.80e-06 -0.2880311 0.0093148 -7.0e-07 -0.0151484 0.0056114 -3.23e-05 0.0106050 NA -0.0083138 0.0009354 12 -315.2432 655.0356 47.133889 5.681881e-11
7 -5.626706 -2.30e-06 -0.1850096 0.0090594 -7.0e-07 -0.0160232 0.0055363 -7.41e-05 0.0081223 -1.6e-06 NA NA 11 -318.9526 660.3691 52.467379 3.947659e-12
13 -5.591850 -4.10e-06 0.0418431 0.0081474 -1.2e-06 -0.0165353 0.0050064 -2.08e-05 NA -3.4e-06 -0.0062482 NA 11 -319.0855 660.6349 52.733197 3.456350e-12
30 -5.630279 8.00e-07 -0.2368065 0.0093399 -9.0e-07 -0.0154785 0.0054981 NA 0.0086443 -1.0e-06 -0.0107548 0.0002452 12 -318.0553 660.6599 52.758122 3.413543e-12
21 -5.596145 -4.60e-06 0.0511733 0.0082013 -1.0e-06 -0.0166104 0.0050502 -4.75e-05 NA -3.3e-06 NA -0.0018999 11 -319.7398 661.9436 54.041824 1.796607e-12
15 -5.623288 -1.20e-06 -0.2342121 0.0090106 -8.0e-07 -0.0156238 0.0055951 -3.32e-05 0.0091760 -1.4e-06 -0.0076945 NA 12 -322.3572 669.2638 61.362018 4.622673e-14
23 -5.625559 -2.30e-06 -0.1802420 0.0090210 -7.0e-07 -0.0159975 0.0055326 -7.14e-05 0.0079155 -1.6e-06 NA -0.0005036 12 -323.1385 670.8264 62.924614 2.116312e-14
29 -5.591361 -4.00e-06 0.0398539 0.0081165 -1.1e-06 -0.0164555 0.0050275 -1.84e-05 NA -3.3e-06 -0.0058413 -0.0010734 12 -323.2685 671.0863 63.184602 1.858335e-14
31 -5.625254 -1.10e-06 -0.2459261 0.0090800 -8.0e-07 -0.0156498 0.0056052 -3.60e-05 0.0096237 -1.4e-06 -0.0081240 0.0009472 13 -326.4928 679.6277 71.725922 2.596570e-16

Run the top model and 100 models for each one:

#In both cases the top model is 1/2/3/4/5/6 no interaction terms. With no models within delta < 4: 

#Run model for DR
Armenta.MCC.DR.top <- gls(MCC.DR ~ Colour.discriminability
                         + log(range.size.m2)
                         + bioclim4 #Seasonal variation
                         + residuals.PC1 #Spatial variation
                         + PC1.LIG #Long-term climate variation
                         + NPP,
                correlation = corPagel(1, phy = pruned.MCC.Armenta.tree, fixed = FALSE), 
                data = MCC.Armenta.model.data, 
                method = "REML")
saveRDS(Armenta.MCC.DR.top, 'data/Armenta.MCC.DR.top.rds')

#Run model for ND
Armenta.MCC.ND.top <- gls(MCC.ND ~ Colour.discriminability
                         + log(range.size.m2)
                         + bioclim4 #Seasonal variation
                         + residuals.PC1 #Spatial variation
                         + PC1.LIG #Long-term climate variation
                         + NPP,
                correlation = corPagel(1, phy = pruned.MCC.Armenta.tree, fixed = TRUE), #lambda = 1
                data = MCC.Armenta.model.data, 
                method = "REML")
saveRDS(Armenta.MCC.ND.top, 'data/Armenta.MCC.ND.top.rds')

#Run the 100 models for DR and ND using the best model:
Armenta.data.noMCC <- inner_join(restricted.data %>% dplyr::select(TipLabel, binomial, SDi, range.size.m2, bioclim4, residuals.PC1, PC1.LIG, NPP), Armenta.data, by = "binomial") %>% filter(Colour.discriminability != "NA")
#Take the restricted data and make it simpler with just responses and predictors.Note that we join the es.values for the 100 trees
Armenta.DR.model.data <- lapply(es.list, function(x) { #es.list is a list of ES values calculated earlier
  left_join(Armenta.data.noMCC %>% dplyr::select(binomial, TipLabel, SDi, Colour.discriminability, range.size.m2, bioclim4, residuals.PC1, PC1.LIG, NPP), 
            x %>% as.data.frame() %>% tibble::rownames_to_column() %>% `colnames<-`(c("TipLabel", "DR")), 
            by = "TipLabel")
})

#PGLS needs tiplabel as rowname
Armenta.DR.model.data <- lapply(Armenta.DR.model.data, function(x) {
  tibble::column_to_rownames(x, "TipLabel")})

#Prune the trees
Armenta.pruned.trees<-lapply(passerine.trees, function(x) {
  drop.tip(x,x$tip.label[-match(MCC.Armenta.model.data$TipLabel, x$tip.label)])
})

#Use mapply to create a list of PGLS global models
Armenta.DR.pgls.models <- mcmapply(function(x,y) {
  gls(DR ~ Colour.discriminability 
         + log(range.size.m2)
         + bioclim4 #Seasonal variation
         + residuals.PC1 #Spatial variation
         + PC1.LIG #Long-term climate variation
         + NPP,
    correlation = corPagel(1, phy = y, fixed = FALSE), 
    data = x, 
    method = "REML")
}, x = Armenta.DR.model.data, y = Armenta.pruned.trees,
SIMPLIFY = FALSE, #Prevents the catastrophic loss of structure, names and attributes
mc.cores = 8) #Specify core number 

saveRDS(Armenta.DR.pgls.models, "data/Armenta.DR.pgls.models.rds")

#Now for Node Density:
Armenta.ND.model.data <- lapply(nd.list, function(x) {
  left_join(Armenta.data.noMCC %>% dplyr::select(binomial, TipLabel, SDi, Colour.discriminability, range.size.m2, bioclim4, residuals.PC1, PC1.LIG, NPP), 
            x %>% as.data.frame() %>% tibble::rownames_to_column() %>% `colnames<-`(c("TipLabel", "ND")), 
            by = "TipLabel")
})

#PGLS needs tiplabel as rowname
Armenta.ND.model.data <- lapply(Armenta.ND.model.data, function(x) {
  tibble::column_to_rownames(x, "TipLabel")})

#Use mapply to create a list of PGLS global models
Armenta.ND.pgls.models <- mcmapply(function(x,y) {
gls(ND ~ Colour.discriminability 
         + log(range.size.m2)
         + bioclim4 #Seasonal variation
         + residuals.PC1 #Spatial variation
         + PC1.LIG #Long-term climate variation
         + NPP,
    corPagel(1, phy = y, fixed = TRUE), 
    data = x, 
    method = "REML")
}, x = Armenta.ND.model.data, y = Armenta.pruned.trees, 
SIMPLIFY = FALSE, #Prevents the catastrophic loss of structure, names and attributes
mc.cores = 8) #Specify core number 

saveRDS(Armenta.ND.pgls.models, "data/Armenta.ND.pgls.models.rds")

#Run the BAMM models

Armenta.MCC.Lambda.top <- gls(mean.lambda ~ Colour.discriminability
                      + log(range.size.m2)
                      + bioclim4 #Seasonal variation
                      + residuals.PC1 #Spatial variation
                      + PC1.LIG #Long-term climate variation
                      + NPP,
                weights = ~ sqrt(var.lambda),
                correlation = corBrownian(phy = pruned.MCC.Armenta.tree), 
                data = MCC.Armenta.model.data, 
                method = "REML")
saveRDS(Armenta.MCC.Lambda.top, 'data/Armenta.MCC.Lambda.top.rds')

Armenta.MCC.Mu.top <- gls(mean.mu ~ Colour.discriminability
                      + log(range.size.m2)
                      + bioclim4 #Seasonal variation
                      + residuals.PC1 #Spatial variation
                      + PC1.LIG #Long-term climate variation
                      + NPP,
                weights = ~ sqrt(var.mu),
                correlation = corBrownian(phy = pruned.MCC.Armenta.tree), 
                data = MCC.Armenta.model.data, 
                method = "REML")
saveRDS(Armenta.MCC.Mu.top, 'data/Armenta.MCC.Mu.top.rds')

Armenta.BAMM.model.data <- lapply(BAMM.df, function(x) { #es.list is a list of ES values calculated earlier
  left_join(Armenta.data.noMCC %>% dplyr::select(binomial, TipLabel, SDi, Colour.discriminability, range.size.m2, bioclim4, residuals.PC1, PC1.LIG, NPP), 
            x %>% as.data.frame(), 
            by = "TipLabel")
})

#PGLS needs tiplabel as rowname
Armenta.BAMM.model.data <- lapply(Armenta.BAMM.model.data, function(x) {
  tibble::column_to_rownames(x, "TipLabel")})

#Use mapply to create a list of PGLS global models
Armenta.BAMM.lambda.pgls.models <- mcmapply(function(x,y) {
  gls(mean.lambda ~ Colour.discriminability 
         + log(range.size.m2)
         + bioclim4 #Seasonal variation
         + residuals.PC1 #Spatial variation
         + PC1.LIG #Long-term climate variation
         + NPP,
    weights = ~ sqrt(var.lambda),      
    corBrownian(phy = y), #lambda = 1
    data = x, 
    method = "REML")
}, x = Armenta.BAMM.model.data, y = Armenta.pruned.trees,
SIMPLIFY = FALSE, #Prevents the catastrophic loss of structure, names and attributes
mc.cores = 8) #Specify core number 

saveRDS(Armenta.BAMM.lambda.pgls.models, "data/Armenta.BAMM.lambda.pgls.models.rds")

#Use mapply to create a list of PGLS global models
Armenta.BAMM.mu.pgls.models <- mcmapply(function(x,y) {
  gls(mean.mu ~ Colour.discriminability
         + log(range.size.m2)
         + bioclim4 #Seasonal variation
         + residuals.PC1 #Spatial variation
         + PC1.LIG #Long-term climate variation
         + NPP,
    weights = ~ sqrt(var.mu),      
    corBrownian(phy = y), #lambda = 1
    data = x, 
    method = "REML")
}, x = Armenta.BAMM.model.data, y = Armenta.pruned.trees,
SIMPLIFY = FALSE, #Prevents the catastrophic loss of structure, names and attributes
mc.cores = 8) #Specify core number 

saveRDS(Armenta.BAMM.mu.pgls.models, "data/Armenta.BAMM.mu.pgls.models.rds")

Using the MCC models as well as the 100 PGLS models we can generate the plots similar to Supp Fig S8:

Armenta.DR.pgls.models <- readRDS("data/Armenta.DR.pgls.models.rds")
Armenta.MCC.DR.top <- readRDS('data/Armenta.MCC.DR.top.rds')

Armenta.DR.pgls.summary <- lapply(Armenta.DR.pgls.models, function(x) {
  data.frame(x$coefficients, confint(x)) %>% tibble::rownames_to_column()
})

Armenta.MCC.DR.summary <- data.frame(Armenta.MCC.DR.top$coefficients, confint(Armenta.MCC.DR.top)) %>% tibble::rownames_to_column()

Armenta.DR.pgls.summary <- bind_rows(Armenta.DR.pgls.summary) %>% `colnames<-`(c("Parameter", "Estimate", "LCI", "UCI"))

colnames(Armenta.MCC.DR.summary) <- c("Parameter", "Estimate", "LCI", "UCI")

parameter_names <- c(
                    `bioclim4` = "Temperature Seasonality",
                    `log(range.size.m2)` = "Range Size (log-transformed)",
                    `NPP` = "NPP",
                    `PC1.LIG` = "Long-term Temperature Variation",
                    `residuals.PC1` = "Spatial Temperature Variation",
                    `Colour.discriminability` = "Sexual Dichromatism"
                    )

Armenta.DR.pgls.summary$Parameter = factor(Armenta.DR.pgls.summary$Parameter, levels=c('bioclim4','log(range.size.m2)','NPP','PC1.LIG', 'residuals.PC1', 'Colour.discriminability'))

Armenta.MCC.DR.summary$Parameter = factor(Armenta.MCC.DR.summary$Parameter, levels=c('bioclim4','log(range.size.m2)','NPP','PC1.LIG', 'residuals.PC1', 'Colour.discriminability'))

Armenta.DR.plot <-Armenta.DR.pgls.summary %>% filter(Parameter != "(Intercept)") %>% 
  ggplot(aes(x = Estimate, y = Parameter, fill = Parameter)) +
  geom_density(aes(y = ..scaled..),
               position = position_nudge(y= 0.6))+
  geom_jitter(shape = 21, alpha = 0.5, size = 0.75)+
  geom_vline(xintercept = 0)+
  theme_minimal()+
  theme(axis.text.y = element_blank(),
        legend.position = "none")

Armenta.DR.plot <- Armenta.DR.plot + geom_errorbarh(data = Armenta.MCC.DR.summary %>% filter(Parameter != "(Intercept)"), 
                 aes(xmin = LCI,
                     xmax = UCI, y = Parameter,
                     color = Parameter), 
                 height = 0, show.legend = F, 
                 position = position_nudge(y= -0.75))+
  
  geom_point(data = Armenta.MCC.DR.summary %>% filter(Parameter != "(Intercept)"), 
             aes(x = Estimate, y = Parameter, fill = Parameter), 
             shape=21, color = "grey20",
             size = 3,
             position = position_nudge(y= -0.75)) +
  scale_y_discrete(expand = c(0.5,.5))+
  facet_wrap(~ Parameter, scales = "free", nrow = 6, labeller = as_labeller(parameter_names))+
  scale_fill_brewer(palette = "Dark2")+
  scale_color_brewer(palette = "Dark2")+
  theme(panel.grid.major.y = element_blank(),
        axis.title.y = element_blank(),
        axis.text.x = element_text(size = 7),
        plot.margin=unit(c(.5,.5,.5,.5),"cm"))+
  ggtitle("Armenta DR")

####################################
#For ND

Armenta.ND.pgls.models <- readRDS("data/Armenta.ND.pgls.models.rds")
Armenta.MCC.ND.top <- readRDS('data/Armenta.MCC.ND.top.rds')

Armenta.ND.pgls.summary <- lapply(Armenta.ND.pgls.models, function(x) {
  data.frame(x$coefficients, confint(x)) %>% tibble::rownames_to_column()
})

Armenta.MCC.ND.summary <- data.frame(Armenta.MCC.ND.top$coefficients, confint(Armenta.MCC.ND.top)) %>% tibble::rownames_to_column()

Armenta.ND.pgls.summary <- bind_rows(Armenta.ND.pgls.summary) %>% `colnames<-`(c("Parameter", "Estimate", "LCI", "UCI"))

colnames(Armenta.MCC.ND.summary) <- c("Parameter", "Estimate", "LCI", "UCI")

Armenta.ND.pgls.summary$Parameter = factor(Armenta.ND.pgls.summary$Parameter, levels=c('bioclim4','log(range.size.m2)','NPP','PC1.LIG', 'residuals.PC1', 'Colour.discriminability'))

Armenta.MCC.ND.summary$Parameter = factor(Armenta.MCC.ND.summary$Parameter, levels=c('bioclim4','log(range.size.m2)','NPP','PC1.LIG', 'residuals.PC1', 'Colour.discriminability'))

Armenta.ND.plot <-Armenta.ND.pgls.summary %>% filter(Parameter != "(Intercept)") %>% 
  ggplot(aes(x = Estimate, y = Parameter, fill = Parameter)) +
  geom_density(aes(y = ..scaled..),
               position = position_nudge(y= 0.6))+
  geom_jitter(shape = 21, alpha = 0.5, size = 0.75)+
  geom_vline(xintercept = 0)+
  theme_minimal()+
  theme(axis.text.y = element_blank(),
        legend.position = "none")

Armenta.ND.plot <- Armenta.ND.plot + geom_errorbarh(data = Armenta.MCC.ND.summary %>% filter(Parameter != "(Intercept)"), 
                 aes(xmin = LCI,
                     xmax = UCI, y = Parameter,
                     color = Parameter), 
                 height = 0, show.legend = F, 
                 position = position_nudge(y= -0.75))+
  
  geom_point(data = Armenta.MCC.ND.summary %>% filter(Parameter != "(Intercept)"), 
             aes(x = Estimate, y = Parameter, fill = Parameter), 
             shape=21, color = "grey20",
             size = 3,
             position = position_nudge(y= -0.75)) +
  scale_y_discrete(expand = c(0.5,.5))+
  facet_wrap(~ Parameter, scales = "free", nrow = 6, labeller = as_labeller(parameter_names))+
  scale_fill_brewer(palette = "Dark2")+
  scale_color_brewer(palette = "Dark2")+
  theme(panel.grid.major.y = element_blank(),
        axis.title.y = element_blank(),
        axis.text.x = element_text(size = 7),
        plot.margin=unit(c(.5,.5,.5,.5),"cm"))+
  ggtitle("Armenta ND")

######################

#For Lambda
Armenta.BAMM.lambda.pgls.models <- readRDS("data/Armenta.BAMM.lambda.pgls.models.rds")
Armenta.MCC.Lambda.top <- readRDS('data/Armenta.MCC.Lambda.top.rds')

Armenta.Lambda.pgls.summary <- lapply(Armenta.BAMM.lambda.pgls.models, function(x) {
  data.frame(x$coefficients, confint(x)) %>% tibble::rownames_to_column()
})

Armenta.MCC.Lambda.summary <- data.frame(Armenta.MCC.Lambda.top$coefficients, confint(Armenta.MCC.Lambda.top)) %>% tibble::rownames_to_column()

Armenta.Lambda.pgls.summary <- bind_rows(Armenta.Lambda.pgls.summary) %>% `colnames<-`(c("Parameter", "Estimate", "LCI", "UCI"))

colnames(Armenta.MCC.Lambda.summary) <- c("Parameter", "Estimate", "LCI", "UCI")

Armenta.Lambda.pgls.summary$Parameter = factor(Armenta.Lambda.pgls.summary$Parameter, levels=c('bioclim4','log(range.size.m2)','NPP','PC1.LIG', 'residuals.PC1', 'Colour.discriminability'))

Armenta.MCC.Lambda.summary$Parameter = factor(Armenta.MCC.Lambda.summary$Parameter, levels=c('bioclim4','log(range.size.m2)','NPP','PC1.LIG', 'residuals.PC1', 'Colour.discriminability'))

Armenta.Lambda.plot <-Armenta.Lambda.pgls.summary %>% filter(Parameter != "(Intercept)") %>% 
  ggplot(aes(x = Estimate, y = Parameter, fill = Parameter)) +
  geom_density(aes(y = ..scaled..),
               position = position_nudge(y= 0.6))+
  geom_jitter(shape = 21, alpha = 0.5, size = 0.75)+
  geom_vline(xintercept = 0)+
  theme_minimal()+
  theme(axis.text.y = element_blank(),
        legend.position = "none")

Armenta.Lambda.plot <- Armenta.Lambda.plot + geom_errorbarh(data = Armenta.MCC.Lambda.summary %>% filter(Parameter != "(Intercept)"), 
                 aes(xmin = LCI,
                     xmax = UCI, y = Parameter,
                     color = Parameter), 
                 height = 0, show.legend = F, 
                 position = position_nudge(y= -0.75))+
  
  geom_point(data = Armenta.MCC.Lambda.summary %>% filter(Parameter != "(Intercept)"), 
             aes(x = Estimate, y = Parameter, fill = Parameter), 
             shape=21, color = "grey20",
             size = 3,
             position = position_nudge(y= -0.75)) +
  scale_y_discrete(expand = c(0.5,.5))+
  facet_wrap(~ Parameter, scales = "free", nrow = 6, labeller = as_labeller(parameter_names))+
  scale_fill_brewer(palette = "Dark2")+
  scale_color_brewer(palette = "Dark2")+
  theme(panel.grid.major.y = element_blank(),
        axis.title.y = element_blank(),
        axis.text.x = element_text(size = 7),
        plot.margin=unit(c(.5,.5,.5,.5),"cm"))+
  ggtitle("Armenta BAMM Speciaton")

######################

#For Mu

Armenta.BAMM.Mu.pgls.models <- readRDS("data/Armenta.BAMM.mu.pgls.models.rds")
Armenta.MCC.Mu.top <- readRDS('data/Armenta.MCC.Mu.top.rds')

Armenta.Mu.pgls.summary <- lapply(Armenta.BAMM.Mu.pgls.models, function(x) {
  data.frame(x$coefficients, confint(x)) %>% tibble::rownames_to_column()
})

Armenta.MCC.Mu.summary <- data.frame(Armenta.MCC.Mu.top$coefficients, confint(Armenta.MCC.Mu.top)) %>% tibble::rownames_to_column()

Armenta.Mu.pgls.summary <- bind_rows(Armenta.Mu.pgls.summary) %>% `colnames<-`(c("Parameter", "Estimate", "LCI", "UCI"))

colnames(Armenta.MCC.Mu.summary) <- c("Parameter", "Estimate", "LCI", "UCI")

Armenta.Mu.pgls.summary$Parameter = factor(Armenta.Mu.pgls.summary$Parameter, levels=c('bioclim4','log(range.size.m2)','NPP','PC1.LIG', 'residuals.PC1', 'Colour.discriminability'))

Armenta.MCC.Mu.summary$Parameter = factor(Armenta.MCC.Mu.summary$Parameter, levels=c('bioclim4','log(range.size.m2)','NPP','PC1.LIG', 'residuals.PC1', 'Colour.discriminability'))

Armenta.Mu.plot <-Armenta.Mu.pgls.summary %>% filter(Parameter != "(Intercept)") %>% 
  ggplot(aes(x = Estimate, y = Parameter, fill = Parameter)) +
  geom_density(aes(y = ..scaled..),
               position = position_nudge(y= 0.6))+
  geom_jitter(shape = 21, alpha = 0.5, size = 0.75)+
  geom_vline(xintercept = 0)+
  theme_minimal()+
  theme(axis.text.y = element_blank(),
        legend.position = "none")

Armenta.Mu.plot <- Armenta.Mu.plot + geom_errorbarh(data = Armenta.MCC.Mu.summary %>% filter(Parameter != "(Intercept)"), 
                 aes(xmin = LCI,
                     xmax = UCI, y = Parameter,
                     color = Parameter), 
                 height = 0, show.legend = F, 
                 position = position_nudge(y= -0.75))+
  
  geom_point(data = Armenta.MCC.Mu.summary %>% filter(Parameter != "(Intercept)"), 
             aes(x = Estimate, y = Parameter, fill = Parameter), 
             shape=21, color = "grey20",
             size = 3,
             position = position_nudge(y= -0.75)) +
  scale_y_discrete(expand = c(0.5,.5))+
  facet_wrap(~ Parameter, scales = "free", nrow = 6, labeller = as_labeller(parameter_names))+
  scale_fill_brewer(palette = "Dark2")+
  scale_color_brewer(palette = "Dark2")+
  theme(panel.grid.major.y = element_blank(),
        axis.title.y = element_blank(),
        axis.text.x = element_text(size = 7),
        plot.margin=unit(c(.5,.5,.5,.5),"cm"))+
  ggtitle("Armenta BAMM Extinction")

Plot the results

grid.arrange(symmetrise_scale(Armenta.DR.plot, "x"),
             symmetrise_scale(Armenta.ND.plot, "x"),
             symmetrise_scale(Armenta.Lambda.plot, "x"), 
             symmetrise_scale(Armenta.Mu.plot, "x"), 
             nrow = 1)

Figure S10: Measures of sexual dichromatism using a more precise measure in spectrophotometry do not change the main patterns seen within the PGLS models analysed here. This dataset is a subset of the complete dataset (n = 581), thus drawing conclusions for the other predictors (e.g. borderline long term temperature variation and spatial temperature variation) potentially risks type I error; additionally these environmental predictors were not measured differently in this analysis compared to the previous above so we have no obvious reason to make conclusions from this subset for the effects of the environmental predictors.

Table S11: MCC estimates from the above plot are presented as numerical values below.

Armenta.MCC.DR.summary %>% filter(Parameter != "(Intercept)") %>% pander(caption = "MCC DR Estimates from Spec measures")
MCC DR Estimates from Spec measures
Parameter Estimate LCI UCI
Colour.discriminability 0.007071 -0.0482 0.06234
log(range.size.m2) 0.009034 -0.02176 0.03982
bioclim4 0.0001789 -8.595e-05 0.0004437
residuals.PC1 0.0156 -0.0073 0.0385
PC1.LIG -0.03157 -0.07051 0.007379
NPP 7.814e-06 -5.345e-06 2.097e-05
Armenta.MCC.ND.summary %>% filter(Parameter != "(Intercept)") %>% pander(caption = "MCC ND Estimates from Spec measures")
MCC ND Estimates from Spec measures
Parameter Estimate LCI UCI
Colour.discriminability 8.285e-05 -0.003516 0.003682
log(range.size.m2) -3.776e-05 -0.001739 0.001663
bioclim4 6.941e-06 -8.315e-06 2.22e-05
residuals.PC1 0.001053 -0.0001951 0.002302
PC1.LIG -0.0008841 -0.003057 0.001288
NPP 4.626e-07 -2.874e-07 1.213e-06
Armenta.MCC.Lambda.summary %>% filter(Parameter != "(Intercept)") %>% pander(caption = "MCC BAMM Speciation Estimates from Spec measures")
MCC BAMM Speciation Estimates from Spec measures
Parameter Estimate LCI UCI
Colour.discriminability -0.0005285 -0.01583 0.01478
log(range.size.m2) 0.001845 -0.005712 0.009402
bioclim4 -2.585e-05 -9.338e-05 4.167e-05
residuals.PC1 0.002489 -0.003265 0.008243
PC1.LIG -0.004244 -0.01356 0.005072
NPP -1.636e-06 -4.926e-06 1.654e-06
Armenta.MCC.Mu.summary %>% filter(Parameter != "(Intercept)") %>% pander(caption = "MCC BAMM Extinction Estimates from Spec measures")
MCC BAMM Extinction Estimates from Spec measures
Parameter Estimate LCI UCI
Colour.discriminability -0.0007063 -0.03742 0.036
log(range.size.m2) 0.009819 -0.007562 0.0272
bioclim4 -2.645e-06 -0.0001577 0.0001524
residuals.PC1 0.004753 -0.008041 0.01755
PC1.LIG -0.01581 -0.03777 0.006146
NPP -6.949e-07 -8.335e-06 6.946e-06

Table S12: Highest posterior density intervals are calculated for the above figure (ie; models using Spectrophotometry values of sexual dichromatism). The HPD range is determined using the hdi function with a 95 % credible interval. These intervals do not take into account the variance associated with each interval and thus are not an estimate of model precision. Intervals not overlapping zero suggest that 95 % of trees from the posterior generate a model estimate for the given parameter that in the same direction (+ or -). These intervals are calculated in the same way as in Table S7.

Armenta.hpd.DR.top <- list()
Armenta.DR.pgls.summary <- na.omit(Armenta.DR.pgls.summary)
for(x in unique(Armenta.DR.pgls.summary$Parameter)) {
Armenta.hpd.DR.top[[x]] = hdi(Armenta.DR.pgls.summary %>% filter(Parameter == x) %>% dplyr::select(Estimate))
}
Armenta.hpd.DR.top <- bind_rows(Armenta.hpd.DR.top) %>% `rownames<-`(c("Lower", "Upper"))


saveRDS(Armenta.hpd.DR.top, 'data/Armenta.hpd.DR.top.rds')

#For ND
Armenta.hpd.ND.top <- list()
Armenta.ND.pgls.summary <- na.omit(Armenta.ND.pgls.summary)
for(x in unique(Armenta.ND.pgls.summary$Parameter)) {
Armenta.hpd.ND.top[[x]] = hdi(Armenta.ND.pgls.summary %>% filter(Parameter == x) %>% dplyr::select(Estimate))
}
Armenta.hpd.ND.top <- bind_rows(Armenta.hpd.ND.top) %>% `rownames<-`(c("Lower", "Upper"))

saveRDS(Armenta.hpd.ND.top, 'data/Armenta.hpd.ND.top.rds')

###############
Armenta.hpd.Lambda.top <- list()
Armenta.Lambda.pgls.summary <- na.omit(Armenta.Lambda.pgls.summary)
for(x in unique(Armenta.Lambda.pgls.summary$Parameter)) {
Armenta.hpd.Lambda.top[[x]] = hdi(Armenta.Lambda.pgls.summary %>% filter(Parameter == x) %>% dplyr::select(Estimate))
}
Armenta.hpd.Lambda.top <- bind_rows(Armenta.hpd.Lambda.top) %>% `rownames<-`(c("Lower", "Upper")) 
saveRDS(Armenta.hpd.Lambda.top, 'data/Armenta.hpd.Lambda.top.rds')

################
Armenta.hpd.Mu.top <- list()
Armenta.Mu.pgls.summary <- na.omit(Armenta.Mu.pgls.summary)
for(x in unique(Armenta.Mu.pgls.summary$Parameter)) {
Armenta.hpd.Mu.top[[x]] = hdi(Armenta.Mu.pgls.summary %>% filter(Parameter == x) %>% dplyr::select(Estimate))
}
Armenta.hpd.Mu.top <- bind_rows(Armenta.hpd.Mu.top) %>% `rownames<-`(c("Lower", "Upper"))

saveRDS(Armenta.hpd.Mu.top, 'data/Armenta.hpd.Mu.top.rds')

Armenta.hpd.DR.top %>% pander(split.table = Inf, digits = 3, caption = "Armenta DR HPD Interval")
Armenta DR HPD Interval
  Colour.discriminability log(range.size.m2) bioclim4 residuals.PC1 PC1.LIG NPP
Lower -0.0178 -0.0148 -0.000139 0.0108 -0.0424 -4.23e-06
Upper 0.0349 0.0316 0.00016 0.0347 0.00792 1.45e-05
Armenta.hpd.ND.top %>% pander(split.table = Inf, digits = 3, caption = "Armenta ND HPD Interval")
Armenta ND HPD Interval
  Colour.discriminability log(range.size.m2) bioclim4 residuals.PC1 PC1.LIG NPP
Lower -0.00176 -0.000596 -4.97e-06 5e-04 -0.0023 -7e-08
Upper 0.00173 0.00133 1.37e-05 0.00184 0.000556 6.84e-07
Armenta.hpd.Lambda.top %>% pander(split.table = Inf, digits = 3, caption = "Armenta BAMM Speciation HPD Interval")
Armenta BAMM Speciation HPD Interval
  Colour.discriminability log(range.size.m2) bioclim4 residuals.PC1 PC1.LIG NPP
Lower -0.0127 -0.00713 -4.41e-05 -0.00289 -0.0173 -4.05e-06
Upper 0.011 0.00837 5.98e-05 0.00677 0.00505 3.69e-06
Armenta.hpd.Mu.top %>% pander(split.table = Inf, digits = 3, caption = "Armenta BAMM Extinction HPD Interval")
Armenta BAMM Extinction HPD Interval
  Colour.discriminability log(range.size.m2) bioclim4 residuals.PC1 PC1.LIG NPP
Lower -0.0282 -0.019 -0.000119 -0.00577 -0.0377 -1.05e-05
Upper 0.0172 0.0207 0.000136 0.0158 0.00783 1.37e-05

Analysis using male-bias measure of sexual selection

Sexual dichromatism is expected to be a good measure of sexual selection strength in birds. However the relationship between sexual dichromatism and male-biased measures of sexual selection (social mating system, sexual size dimorphism, and paternal care) is expected to be relatively noisy given the precision of the measurment used (Dale et al. 2015). Here we use a dataset of sexual selection for 2,465 species within the orginal dataset of sexual dichromatism scores from RGB measures. This male-bias sexual selection score is based on three components, combined in a phylogenetic pca (ppca) with the following loadings:

  • sexual size dimorphism: 0.37
  • social polygyny: 0.57
  • paternal care: -0.57

As such, high values for this score are expected to have high sexual selection (high dimorphism, high polygyny and low paternal care). PGLS models using this dataset are calculated in the same process as above.

SS.subset <- plumage.scores %>% filter(Sexual_selection_ppca != "NA")
SS.subset <- inner_join(SS.subset, restricted.data %>% dplyr::select(TipLabel, binomial, range.size.m2, bioclim4, residuals.PC1, PC1.LIG, NPP, MCC.DR, MCC.ND), by = "TipLabel")
SS.subset <- inner_join(SS.subset, MCC.BAMM.df, by = "TipLabel")

SS.subset %>% ggplot(aes(x = SDi, y = Sexual_selection_ppca))+
  geom_point()+
  geom_smooth(method = "lm")+
  theme_minimal()

Figure S11: The relationship between absolute sexual dichromatism and the sexual selection. Although positively correlated the distribution is very noisy with species having high sexual dichromatism and low sexual selection and visa versa.Colour.discriminability

data_frame(
R2 = summary(lm(SDi ~ Sexual_selection_ppca,
   data = SS.subset))$r.squared,
r = cor(SS.subset$SDi, SS.subset$Sexual_selection_ppca)) %>% pander(digits = 2)
R2 r
0.12 0.34
#Prune tree
pruned.MCC.Subset.tree <- drop.tip(MCC.passerine,MCC.passerine$tip.label[-match(SS.subset$TipLabel, MCC.passerine$tip.label)])

saveRDS(pruned.MCC.Subset.tree, 'data/pruned.MCC.Subset.tree.rds')

#Set rownames to match tree
rownames(SS.subset) <- SS.subset$TipLabel

#Run a corPagel model to estimate lambda for DR
MCC.DR.Subset.global <- gls(MCC.DR ~ Sexual_selection_ppca
                         + log(range.size.m2)
                         + bioclim4 #Seasonal variation
                         + residuals.PC1 #Spatial variation
                         + PC1.LIG #Long-term climate variation
                         + NPP
                         + Sexual_selection_ppca*log(range.size.m2)
                         + Sexual_selection_ppca*bioclim4
                         + Sexual_selection_ppca*residuals.PC1
                         + Sexual_selection_ppca*PC1.LIG
                         + Sexual_selection_ppca*NPP,
                correlation = corPagel(1, phy = pruned.MCC.Subset.tree, fixed = FALSE), 
                data = SS.subset, 
                method = "REML")

#Run a corPagel model to estimate lambda for DR
MCC.ND.Subset.global <- gls(MCC.ND ~ Sexual_selection_ppca
                         + log(range.size.m2)
                         + bioclim4 #Seasonal variation
                         + residuals.PC1 #Spatial variation
                         + PC1.LIG #Long-term climate variation
                         + NPP
                         + Sexual_selection_ppca*log(range.size.m2)
                         + Sexual_selection_ppca*bioclim4
                         + Sexual_selection_ppca*residuals.PC1
                         + Sexual_selection_ppca*PC1.LIG
                         + Sexual_selection_ppca*NPP,
                correlation = corPagel(1, phy = pruned.MCC.Subset.tree, fixed = TRUE), #lambda = 1
                data = SS.subset, 
                method = "REML")

#Run a corPagel model to estimate lambda for DR
MCC.Lambda.Subset.global <- gls(mean.lambda ~ Sexual_selection_ppca
                         + log(range.size.m2)
                         + bioclim4 #Seasonal variation
                         + residuals.PC1 #Spatial variation
                         + PC1.LIG #Long-term climate variatioColour.discriminability
                         + NPP
                         + Sexual_selection_ppca*log(range.size.m2)
                         + Sexual_selection_ppca*bioclim4
                         + Sexual_selection_ppca*residuals.PC1
                         + Sexual_selection_ppca*PC1.LIG
                         + Sexual_selection_ppca*NPP,
                weights = ~ sqrt(var.lambda),
                correlation = corPagel(1, phy = pruned.MCC.Subset.tree, fixed = TRUE), 
                data = SS.subset, 
                method = "REML")

#Run a corPagel model to estimate lambda for DR
MCC.Extinction.Subset.global <- gls(mean.mu ~ Sexual_selection_ppca
                         + log(range.size.m2)
                         + bioclim4 #Seasonal variation
                         + residuals.PC1 #Spatial variation
                         + PC1.LIG #Long-term climate variation
                         + NPP
                         + Sexual_selection_ppca*log(range.size.m2)
                         + Sexual_selection_ppca*bioclim4
                         + Sexual_selection_ppca*residuals.PC1
                         + Sexual_selection_ppca*PC1.LIG
                         + Sexual_selection_ppca*NPP,
                weights = ~ sqrt(var.mu),
                correlation = corPagel(1, phy = pruned.MCC.Subset.tree, fixed = TRUE), 
                data = SS.subset, 
                method = "REML")
#Set up cluster
cores<-8
clusterType <- if(length(find.package("snow", quiet = TRUE))) "SOCK" else "PSOCK"
clust <- try(makeCluster(getOption("cl.cores", cores), type = clusterType))
myclust<-clust

#Export data and packages to cluster
clusterExport(myclust, c("SS.subset"), envir=environment())
clusterExport(myclust, c("pruned.MCC.Subset.tree"), envir=environment())
clusterEvalQ(myclust, library(nlme))
clusterEvalQ(myclust, library(ape))
clusterEvalQ(myclust, library(MuMIn))

#Dredged models:

Subset.dredged.ND.model <- pdredge(MCC.ND.Subset.global, fixed = c("Sexual_selection_ppca", "log(range.size.m2)", "bioclim4", "residuals.PC1", "PC1.LIG", "NPP"), cluster=myclust)

saveRDS(Subset.dredged.ND.model, 'data/Subset.dredged.ND.model.rds')

Subset.dredged.DR.model <- pdredge(MCC.DR.Subset.global, fixed = c("Sexual_selection_ppca", "log(range.size.m2)", "bioclim4", "residuals.PC1", "PC1.LIG", "NPP"), cluster=myclust)

saveRDS(Subset.dredged.DR.model, 'data/Subset.dredged.DR.model.rds')

Subset.dredged.spec.model <- pdredge(MCC.Lambda.Subset.global, fixed = c("Sexual_selection_ppca", "log(range.size.m2)", "bioclim4", "residuals.PC1", "PC1.LIG", "NPP"), cluster=myclust)

saveRDS(Subset.dredged.spec.model, 'data/Subset.dredged.spec.model.rds')

Subset.dredged.extinct.model <- pdredge(MCC.Extinction.Subset.global, fixed = c("Sexual_selection_ppca", "log(range.size.m2)", "bioclim4", "residuals.PC1", "PC1.LIG", "NPP"), cluster=myclust)

saveRDS(Subset.dredged.spec.model, 'data/Subset.dredged.extinct.model.rds')

Table S13: All top models (\(\delta AICc > 4\)) using male-biased sexual selection measures do not contain interactions between sexual selection and the environmental variables. Thus interaction terms are not included in further analysis.

Subset.dredged.ND.model <- readRDS("data/Subset.dredged.ND.model.rds")
Subset.dredged.DR.model <- readRDS("data/Subset.dredged.DR.model.rds")
Subset.dredged.spec.model <- readRDS("data/Subset.dredged.spec.model.rds")
Subset.dredged.extinct.model <- readRDS("data/Subset.dredged.extinct.model.rds")


kable(Subset.dredged.ND.model, "html", caption = "ND Male-bias Sexual Selection Dredge Table") %>%
  kable_styling() %>%
  scroll_box(width = "100%", height = "500px")
ND Male-bias Sexual Selection Dredge Table
(Intercept) bioclim4 log(range.size.m2) NPP PC1.LIG residuals.PC1 Sexual_selection_ppca bioclim4:Sexual_selection_ppca log(range.size.m2):Sexual_selection_ppca NPP:Sexual_selection_ppca PC1.LIG:Sexual_selection_ppca residuals.PC1:Sexual_selection_ppca df logLik AICc delta weight
0 0.1484807 -5e-07 -0.0001800 0 0.0001073 0.0000788 0.0004381 NA NA NA NA NA 8 5451.548 -10887.04 0.00000 9.994587e-01
16 0.1482246 -5e-07 -0.0001720 0 0.0000929 0.0000644 0.0003782 NA NA NA NA -0.0001167 9 5444.064 -10870.05 16.98341 2.050523e-04
8 0.1479484 -5e-07 -0.0001570 0 0.0000777 0.0001046 0.0004928 NA NA NA 0.0001280 NA 9 5444.026 -10869.98 17.05889 1.974574e-04
2 0.1483243 -5e-07 -0.0001765 0 0.0001017 0.0000834 0.0018018 NA -4.87e-05 NA NA NA 9 5443.645 -10869.22 17.82105 1.348878e-04
1 0.1478764 -9e-07 -0.0001561 0 0.0000738 0.0000871 0.0012068 -2.3e-06 NA NA NA NA 9 5440.065 -10862.06 24.98016 3.761776e-06
4 0.1484220 -6e-07 -0.0001730 0 0.0001024 0.0000814 0.0008188 NA NA 0e+00 NA NA 9 5436.296 -10854.52 32.51872 8.677880e-08
24 0.1477981 -5e-07 -0.0001528 0 0.0000692 0.0000925 0.0004458 NA NA NA 0.0001196 -0.0000845 10 5436.455 -10852.82 34.21679 3.712640e-08
18 0.1479767 -6e-07 -0.0001662 0 0.0000832 0.0000688 0.0022340 NA -6.65e-05 NA NA -0.0001323 10 5436.226 -10852.36 34.67399 2.953948e-08
10 0.1478037 -5e-07 -0.0001538 0 0.0000725 0.0001088 0.0017783 NA -4.59e-05 NA 0.0001274 NA 10 5436.116 -10852.14 34.89512 2.644759e-08
17 0.1478092 -9e-07 -0.0001543 0 0.0000700 0.0000806 0.0011337 -2.2e-06 NA NA NA -0.0000480 10 5432.459 -10844.83 42.20969 6.824126e-10
9 0.1474813 -8e-07 -0.0001389 0 0.0000518 0.0001085 0.0011728 -2.1e-06 NA NA 0.0001104 NA 10 5432.373 -10844.66 42.38150 6.262388e-10
3 0.1478909 -9e-07 -0.0001563 0 0.0000742 0.0000866 0.0010291 -2.3e-06 6.70e-06 NA NA NA 10 5432.150 -10844.21 42.82661 5.012840e-10
20 0.1482076 -6e-07 -0.0001673 0 0.0000906 0.0000682 0.0006967 NA NA 0e+00 NA -0.0001026 10 5428.776 -10837.46 49.57473 1.716915e-11
12 0.1478969 -6e-07 -0.0001505 0 0.0000733 0.0001068 0.0008569 NA NA 0e+00 0.0001269 NA 10 5428.762 -10837.43 49.60299 1.692825e-11
6 0.1482266 -6e-07 -0.0001680 0 0.0000951 0.0000872 0.0025093 NA -5.90e-05 0e+00 NA NA 10 5428.428 -10836.77 50.27032 1.212563e-11
26 0.1475846 -6e-07 -0.0001479 0 0.0000610 0.0000958 0.0021036 NA -5.95e-05 NA 0.0001173 -0.0000990 11 5428.596 -10835.09 51.95195 5.230492e-12
5 0.1477933 -1e-06 -0.0001475 0 0.0000674 0.0000903 0.0016611 -2.4e-06 NA -1e-07 NA NA 10 5424.857 -10829.62 57.41215 3.410950e-13
25 0.1474541 -8e-07 -0.0001383 0 0.0000503 0.0001049 0.0011368 -2.0e-06 NA NA 0.0001086 -0.0000240 11 5424.751 -10827.39 59.64290 1.118085e-13
19 0.1477962 -9e-07 -0.0001541 0 0.0000696 0.0000808 0.0012568 -2.1e-06 -4.70e-06 NA NA -0.0000499 11 5424.591 -10827.07 59.96353 9.524668e-14
11 0.1474880 -8e-07 -0.0001390 0 0.0000520 0.0001082 0.0010942 -2.1e-06 2.90e-06 NA 0.0001103 NA 11 5424.457 -10826.81 60.23059 8.334079e-14
28 0.1477791 -6e-07 -0.0001479 0 0.0000668 0.0000965 0.0007717 NA NA 0e+00 0.0001201 -0.0000699 11 5421.172 -10820.24 66.80118 3.119414e-15
22 0.1479328 -7e-07 -0.0001603 0 0.0000797 0.0000735 0.0027753 NA -7.32e-05 0e+00 NA -0.0001181 11 5420.966 -10819.82 67.21289 2.539046e-15
14 0.1477158 -6e-07 -0.0001460 0 0.0000665 0.0001122 0.0024557 NA -5.58e-05 0e+00 0.0001260 NA 11 5420.885 -10819.66 67.37476 2.341650e-15
21 0.1477617 -1e-06 -0.0001470 0 0.0000656 0.0000868 0.0016033 -2.3e-06 NA 0e+00 NA -0.0000251 11 5417.249 -10812.39 74.64696 6.171400e-17
13 0.1474084 -1e-06 -0.0001309 0 0.0000460 0.0001112 0.0016076 -2.1e-06 NA 0e+00 0.0001085 NA 11 5417.149 -10812.19 74.84702 5.583943e-17
7 0.1477851 -1e-06 -0.0001473 0 0.0000671 0.0000906 0.0017584 -2.4e-06 -3.60e-06 -1e-07 NA NA 11 5416.948 -10811.79 75.24841 4.568560e-17
27 0.1474468 -8e-07 -0.0001381 0 0.0000501 0.0001050 0.0012077 -2.0e-06 -2.70e-06 NA 0.0001085 -0.0000251 12 5416.883 -10809.64 77.39880 1.558939e-17
30 0.1475395 -6e-07 -0.0001419 0 0.0000575 0.0001007 0.0026498 NA -6.62e-05 0e+00 0.0001176 -0.0000846 12 5413.338 -10802.55 84.48797 4.502300e-19
29 0.1474072 -1e-06 -0.0001309 0 0.0000459 0.0001110 0.0016049 -2.1e-06 NA 0e+00 0.0001084 -0.0000012 12 5409.540 -10794.95 92.08377 1.009318e-20
23 0.1477348 -1e-06 -0.0001465 0 0.0000647 0.0000871 0.0018596 -2.2e-06 -9.70e-06 0e+00 NA -0.0000288 12 5409.384 -10794.64 92.39589 8.634816e-21
15 0.1473924 -1e-06 -0.0001306 0 0.0000455 0.0001118 0.0017935 -2.1e-06 -6.80e-06 0e+00 0.0001087 NA 12 5409.241 -10794.35 92.68281 7.480788e-21
31 0.1473862 -1e-06 -0.0001305 0 0.0000452 0.0001112 0.0018080 -2.1e-06 -7.70e-06 0e+00 0.0001084 -0.0000042 13 5401.675 -10777.20 109.83562 1.410139e-24
kable(Subset.dredged.DR.model, "html", caption = "DR Male-bias Sexual Selection Dredge Table") %>%
  kable_styling() %>%
  scroll_box(width = "100%", height = "500px")
DR Male-bias Sexual Selection Dredge Table
(Intercept) bioclim4 log(range.size.m2) NPP PC1.LIG residuals.PC1 Sexual_selection_ppca bioclim4:Sexual_selection_ppca log(range.size.m2):Sexual_selection_ppca NPP:Sexual_selection_ppca PC1.LIG:Sexual_selection_ppca residuals.PC1:Sexual_selection_ppca df logLik AICc delta weight
0 -2.644369 2.45e-05 -0.0113055 -7.0e-07 0.0068813 0.0004409 0.0388722 NA NA NA NA NA 9 -1717.437 3452.947 0.00000 9.851357e-01
16 -2.639661 2.53e-05 -0.0115625 -7.0e-07 0.0066554 -0.0000035 0.0381221 NA NA NA NA -0.0043664 10 -1721.450 3462.989 10.04197 6.499963e-03
2 -2.651446 2.27e-05 -0.0111253 -7.0e-07 0.0068731 0.0003113 0.1032726 NA -0.0023191 NA NA NA 10 -1721.901 3463.892 10.94493 4.138415e-03
8 -2.643685 2.47e-05 -0.0113672 -8.0e-07 0.0067933 0.0002826 0.0391218 NA NA NA -0.0017643 NA 10 -1721.915 3463.919 10.97215 4.082471e-03
1 -2.644777 1.73e-05 -0.0112474 -8.0e-07 0.0069016 0.0002884 0.0537640 -4.40e-05 NA NA NA NA 10 -1726.048 3472.186 19.23901 6.543319e-05
18 -2.648357 2.30e-05 -0.0113613 -7.0e-07 0.0066118 -0.0002393 0.1227636 NA -0.0030524 NA NA -0.0050116 11 -1725.773 3473.654 20.70712 3.140527e-05
24 -2.639316 2.54e-05 -0.0115972 -7.0e-07 0.0066022 -0.0000984 0.0383393 NA NA NA -0.0012567 -0.0041819 11 -1725.963 3474.034 21.08668 2.597662e-05
10 -2.650479 2.30e-05 -0.0111903 -7.0e-07 0.0067936 0.0001735 0.1006675 NA -0.0022170 NA -0.0015983 NA 11 -1726.392 3474.892 21.94541 1.690876e-05
4 -2.639136 2.37e-05 -0.0113014 -1.1e-06 0.0067364 0.0005491 0.0550570 NA NA -1.6e-06 NA NA 10 -1729.238 3478.565 25.61792 2.695472e-06
17 -2.641237 1.95e-05 -0.0114476 -8.0e-07 0.0067343 0.0000040 0.0497120 -3.36e-05 NA NA NA -0.0031418 11 -1730.243 3482.593 29.64645 3.596263e-07
9 -2.644302 1.77e-05 -0.0112841 -8.0e-07 0.0068566 0.0002153 0.0532073 -4.20e-05 NA NA -0.0008810 NA 11 -1730.569 3483.246 30.29926 2.594752e-07
3 -2.648526 1.70e-05 -0.0111570 -7.0e-07 0.0068950 0.0002315 0.0876844 -3.97e-05 -0.0012746 NA NA NA 11 -1730.598 3483.303 30.35634 2.521738e-07
26 -2.647793 2.32e-05 -0.0113950 -7.0e-07 0.0065727 -0.0003052 0.1206409 NA -0.0029697 NA -0.0009490 -0.0048550 12 -1730.301 3484.728 31.78144 1.236643e-07
20 -2.635570 2.45e-05 -0.0115406 -1.0e-06 0.0065504 0.0001210 0.0518306 NA NA -1.4e-06 NA -0.0040359 11 -1733.315 3488.737 35.78960 1.666795e-08
6 -2.646854 2.14e-05 -0.0110830 -1.1e-06 0.0067038 0.0004113 0.1349625 NA -0.0027849 -1.9e-06 NA NA 11 -1733.615 3489.338 36.39048 1.234251e-08
12 -2.638336 2.39e-05 -0.0113656 -1.1e-06 0.0066414 0.0003852 0.0556198 NA NA -1.7e-06 -0.0018515 NA 11 -1733.706 3489.519 36.57214 1.127085e-08
19 -2.647224 1.97e-05 -0.0113350 -7.0e-07 0.0066756 -0.0001778 0.1100270 -2.27e-05 -0.0023104 NA NA -0.0040281 12 -1734.642 3493.410 40.46326 1.610686e-09
25 -2.640864 1.99e-05 -0.0114754 -8.0e-07 0.0067004 -0.0000514 0.0493311 -3.21e-05 NA NA -0.0007139 -0.0030923 12 -1734.769 3493.665 40.71787 1.418151e-09
11 -2.648042 1.75e-05 -0.0111939 -8.0e-07 0.0068500 0.0001586 0.0870861 -3.76e-05 -0.0012730 NA -0.0008797 NA 12 -1735.118 3494.364 41.41693 9.998225e-10
5 -2.639640 1.65e-05 -0.0112442 -1.1e-06 0.0067596 0.0003957 0.0694708 -4.36e-05 NA -1.6e-06 NA NA 11 -1737.859 3497.825 44.87789 1.771687e-10
22 -2.644516 2.18e-05 -0.0113100 -1.1e-06 0.0064809 -0.0001171 0.1492103 NA -0.0034126 -1.7e-06 NA -0.0046899 12 -1737.559 3499.245 46.29775 8.710986e-11
28 -2.635111 2.47e-05 -0.0115772 -1.1e-06 0.0064896 0.0000208 0.0524230 NA NA -1.4e-06 -0.0013747 -0.0038249 12 -1737.818 3499.763 46.81566 6.723639e-11
14 -2.645817 2.17e-05 -0.0111498 -1.2e-06 0.0066197 0.0002691 0.1325514 NA -0.0026832 -1.9e-06 -0.0016616 NA 12 -1738.100 3500.327 47.37978 5.071161e-11
27 -2.646817 2.00e-05 -0.0113622 -7.0e-07 0.0066447 -0.0002281 0.1093216 -2.14e-05 -0.0022967 NA -0.0006587 -0.0039771 13 -1739.168 3504.485 51.53820 6.340417e-12
21 -2.637030 1.86e-05 -0.0114218 -1.1e-06 0.0066273 0.0001336 0.0642539 -3.45e-05 NA -1.4e-06 NA -0.0027665 12 -1742.094 3508.315 55.36799 9.343053e-13
7 -2.644376 1.61e-05 -0.0111162 -1.1e-06 0.0067357 0.0003290 0.1181574 -3.75e-05 -0.0017672 -1.8e-06 NA NA 12 -1742.351 3508.828 55.88132 7.228072e-13
13 -2.639083 1.70e-05 -0.0112841 -1.1e-06 0.0067082 0.0003160 0.0690252 -4.13e-05 NA -1.6e-06 -0.0009791 NA 12 -1742.374 3508.875 55.92791 7.061633e-13
30 -2.643870 2.20e-05 -0.0113459 -1.1e-06 0.0064361 -0.0001880 0.1471758 NA -0.0033252 -1.7e-06 -0.0010505 -0.0045120 13 -1742.080 3510.308 57.36050 3.450013e-13
23 -2.643470 1.87e-05 -0.0112857 -1.1e-06 0.0065443 -0.0000596 0.1365454 -2.19e-05 -0.0026920 -1.6e-06 NA -0.0037463 13 -1746.434 3519.016 66.06877 4.434500e-15
29 -2.636582 1.90e-05 -0.0114524 -1.1e-06 0.0065870 0.0000719 0.0640021 -3.28e-05 NA -1.4e-06 -0.0008236 -0.0027038 13 -1746.614 3519.376 66.42937 3.702894e-15
15 -2.643820 1.66e-05 -0.0111564 -1.1e-06 0.0066838 0.0002487 0.1177972 -3.52e-05 -0.0017704 -1.8e-06 -0.0009854 NA 13 -1746.865 3519.878 66.93090 2.881620e-15
31 -2.642995 1.90e-05 -0.0113159 -1.1e-06 0.0065069 -0.0001168 0.1359737 -2.04e-05 -0.0026794 -1.6e-06 -0.0007727 -0.0036828 14 -1750.955 3530.081 77.13422 1.753934e-17
kable(Subset.dredged.spec.model, "html", caption = "BAMM Speciation Male-bias Sexual Selection Dredge Table") %>%
  kable_styling() %>%
  scroll_box(width = "100%", height = "500px")
BAMM Speciation Male-bias Sexual Selection Dredge Table
(Intercept) bioclim4 log(range.size.m2) NPP PC1.LIG residuals.PC1 Sexual_selection_ppca bioclim4:Sexual_selection_ppca log(range.size.m2):Sexual_selection_ppca NPP:Sexual_selection_ppca PC1.LIG:Sexual_selection_ppca residuals.PC1:Sexual_selection_ppca df logLik AICc delta weight
0 -2.291100 -7.6e-06 0.0005471 -1e-07 -0.0005670 0.0004011 0.0009422 NA NA NA NA NA 8 2239.623 -4463.188 0.00000 9.982421e-01
16 -2.291695 -7.7e-06 0.0005655 -1e-07 -0.0005911 0.0003593 0.0008204 NA NA NA NA -0.0002658 9 2233.363 -4448.652 14.53554 6.964394e-04
8 -2.292345 -7.7e-06 0.0005965 -1e-07 -0.0006020 0.0004437 0.0009637 NA NA NA 0.0003152 NA 9 2233.113 -4448.152 15.03575 5.423296e-04
2 -2.290976 -7.6e-06 0.0005418 -1e-07 -0.0005524 0.0003885 -0.0024843 NA 0.0001231 NA NA NA 9 2233.054 -4448.035 15.15316 5.114101e-04
1 -2.291399 -7.9e-06 0.0005605 -1e-07 -0.0005809 0.0004077 0.0013441 -1.2e-06 NA NA NA NA 9 2228.666 -4439.258 23.92997 6.351986e-06
24 -2.292817 -7.8e-06 0.0006107 -1e-07 -0.0006217 0.0004057 0.0008568 NA NA NA 0.0003034 -0.0002315 10 2226.838 -4433.586 29.60133 3.727258e-07
18 -2.291559 -7.6e-06 0.0005602 -1e-07 -0.0005786 0.0003533 -0.0016797 NA 0.0000902 NA NA -0.0002452 10 2226.797 -4433.504 29.68333 3.577527e-07
4 -2.291653 -7.9e-06 0.0005740 -1e-07 -0.0005907 0.0004065 0.0023362 NA NA -2e-07 NA NA 9 2225.687 -4433.301 29.88723 3.230772e-07
10 -2.292220 -7.7e-06 0.0005912 -1e-07 -0.0005877 0.0004314 -0.0023393 NA 0.0001187 NA 0.0003140 NA 10 2226.542 -4432.994 30.19380 2.771635e-07
17 -2.291756 -7.8e-06 0.0005687 -1e-07 -0.0005941 0.0003635 0.0009507 -4.0e-07 NA NA NA -0.0002518 10 2222.457 -4424.824 38.36346 4.663555e-09
3 -2.291393 -7.9e-06 0.0005605 -1e-07 -0.0005682 0.0003937 -0.0032200 -1.9e-06 0.0001720 NA NA NA 10 2222.169 -4424.249 38.93904 3.497280e-09
9 -2.292576 -7.9e-06 0.0006068 -1e-07 -0.0006129 0.0004487 0.0012922 -1.0e-06 NA NA 0.0003117 NA 10 2222.151 -4424.212 38.97548 3.434136e-09
20 -2.292081 -7.9e-06 0.0005863 -1e-07 -0.0006079 0.0003718 0.0020931 NA NA -1e-07 NA -0.0002167 10 2219.415 -4418.741 44.44700 2.226858e-10
26 -2.292680 -7.7e-06 0.0006054 -1e-07 -0.0006090 0.0003997 -0.0016523 NA 0.0000905 NA 0.0003035 -0.0002108 11 2220.273 -4418.437 44.75033 1.913486e-10
12 -2.292870 -8.0e-06 0.0006222 -1e-07 -0.0006246 0.0004485 0.0023252 NA NA -1e-07 0.0003115 NA 10 2219.171 -4418.253 44.93462 1.745048e-10
6 -2.291537 -7.8e-06 0.0005690 -1e-07 -0.0005789 0.0003970 -0.0002205 NA 0.0000895 -1e-07 NA NA 10 2219.115 -4418.140 45.04770 1.649120e-10
19 -2.291677 -7.8e-06 0.0005670 -1e-07 -0.0005823 0.0003628 -0.0022333 -1.1e-06 0.0001231 NA NA -0.0001991 11 2215.992 -4409.877 53.31107 2.647794e-12
25 -2.292858 -7.8e-06 0.0006128 -1e-07 -0.0006237 0.0004086 0.0009470 -3.0e-07 NA NA 0.0003030 -0.0002218 11 2215.932 -4409.756 53.43142 2.493164e-12
5 -2.292014 -8.2e-06 0.0005902 -1e-07 -0.0006074 0.0004143 0.0028381 -1.4e-06 NA -2e-07 NA NA 10 2214.737 -4409.384 53.80399 2.069413e-12
11 -2.292555 -8.0e-06 0.0006063 -1e-07 -0.0006006 0.0004350 -0.0029759 -1.7e-06 0.0001609 NA 0.0003077 NA 11 2215.650 -4409.191 53.99633 1.879667e-12
28 -2.293197 -8.0e-06 0.0006311 -1e-07 -0.0006382 0.0004180 0.0021203 NA NA -1e-07 0.0003025 -0.0001829 11 2212.889 -4403.671 59.51672 1.189444e-13
22 -2.291971 -7.9e-06 0.0005818 -1e-07 -0.0005983 0.0003671 0.0002469 NA 0.0000652 -1e-07 NA -0.0002034 11 2212.849 -4403.591 59.59667 1.142835e-13
14 -2.292756 -7.9e-06 0.0006172 -1e-07 -0.0006133 0.0004393 -0.0001273 NA 0.0000859 -1e-07 0.0003107 NA 11 2212.598 -4403.089 60.09851 8.892208e-14
21 -2.292218 -8.1e-06 0.0005934 -1e-07 -0.0006146 0.0003809 0.0024002 -8.0e-07 NA -1e-07 NA -0.0001865 11 2208.522 -4394.936 68.25187 1.508452e-15
27 -2.292778 -7.9e-06 0.0006110 -1e-07 -0.0006121 0.0004077 -0.0021330 -9.0e-07 0.0001191 NA 0.0003020 -0.0001710 12 2209.466 -4394.805 68.38296 1.412746e-15
7 -2.291971 -8.2e-06 0.0005884 -1e-07 -0.0005954 0.0004025 -0.0009542 -2.0e-06 0.0001395 -1e-07 NA NA 11 2208.232 -4394.357 68.83056 1.129455e-15
13 -2.293157 -8.2e-06 0.0006351 -1e-07 -0.0006382 0.0004545 0.0027477 -1.2e-06 NA -2e-07 0.0003072 NA 11 2208.216 -4394.325 68.86300 1.111284e-15
30 -2.293086 -8.0e-06 0.0006266 -1e-07 -0.0006285 0.0004133 0.0002597 NA 0.0000657 -1e-07 0.0003026 -0.0001694 12 2206.324 -4388.520 74.66758 6.100666e-17
23 -2.292138 -8.1e-06 0.0005912 -1e-07 -0.0006038 0.0003798 -0.0003934 -1.4e-06 0.0001065 -1e-07 NA -0.0001426 12 2202.054 -4379.981 83.20700 8.532288e-19
29 -2.293309 -8.1e-06 0.0006369 -1e-07 -0.0006437 0.0004255 0.0023796 -7.0e-07 NA -1e-07 0.0003014 -0.0001575 12 2201.995 -4379.862 83.32553 8.041324e-19
15 -2.293107 -8.3e-06 0.0006330 -1e-07 -0.0006269 0.0004432 -0.0007647 -1.7e-06 0.0001292 -1e-07 0.0003042 NA 12 2201.708 -4379.289 83.89845 6.038371e-19
31 -2.293228 -8.2e-06 0.0006348 -1e-07 -0.0006332 0.0004243 -0.0003138 -1.2e-06 0.0001027 -1e-07 0.0003006 -0.0001152 13 2195.526 -4364.904 98.28412 4.540571e-22
kable(Subset.dredged.extinct.model, "html", caption = "BAMM Extinction Male-bias Sexual Selection Dredge Table") %>%
  kable_styling() %>%
  scroll_box(width = "100%", height = "500px")
BAMM Extinction Male-bias Sexual Selection Dredge Table
(Intercept) bioclim4 log(range.size.m2) NPP PC1.LIG residuals.PC1 Sexual_selection_ppca bioclim4:Sexual_selection_ppca log(range.size.m2):Sexual_selection_ppca NPP:Sexual_selection_ppca PC1.LIG:Sexual_selection_ppca residuals.PC1:Sexual_selection_ppca df logLik AICc delta weight
0 -2.291100 -7.6e-06 0.0005471 -1e-07 -0.0005670 0.0004011 0.0009422 NA NA NA NA NA 8 2239.623 -4463.188 0.00000 9.982421e-01
16 -2.291695 -7.7e-06 0.0005655 -1e-07 -0.0005911 0.0003593 0.0008204 NA NA NA NA -0.0002658 9 2233.363 -4448.652 14.53554 6.964394e-04
8 -2.292345 -7.7e-06 0.0005965 -1e-07 -0.0006020 0.0004437 0.0009637 NA NA NA 0.0003152 NA 9 2233.113 -4448.152 15.03575 5.423296e-04
2 -2.290976 -7.6e-06 0.0005418 -1e-07 -0.0005524 0.0003885 -0.0024843 NA 0.0001231 NA NA NA 9 2233.054 -4448.035 15.15316 5.114101e-04
1 -2.291399 -7.9e-06 0.0005605 -1e-07 -0.0005809 0.0004077 0.0013441 -1.2e-06 NA NA NA NA 9 2228.666 -4439.258 23.92997 6.351986e-06
24 -2.292817 -7.8e-06 0.0006107 -1e-07 -0.0006217 0.0004057 0.0008568 NA NA NA 0.0003034 -0.0002315 10 2226.838 -4433.586 29.60133 3.727258e-07
18 -2.291559 -7.6e-06 0.0005602 -1e-07 -0.0005786 0.0003533 -0.0016797 NA 0.0000902 NA NA -0.0002452 10 2226.797 -4433.504 29.68333 3.577527e-07
4 -2.291653 -7.9e-06 0.0005740 -1e-07 -0.0005907 0.0004065 0.0023362 NA NA -2e-07 NA NA 9 2225.687 -4433.301 29.88723 3.230772e-07
10 -2.292220 -7.7e-06 0.0005912 -1e-07 -0.0005877 0.0004314 -0.0023393 NA 0.0001187 NA 0.0003140 NA 10 2226.542 -4432.994 30.19380 2.771635e-07
17 -2.291756 -7.8e-06 0.0005687 -1e-07 -0.0005941 0.0003635 0.0009507 -4.0e-07 NA NA NA -0.0002518 10 2222.457 -4424.824 38.36346 4.663555e-09
3 -2.291393 -7.9e-06 0.0005605 -1e-07 -0.0005682 0.0003937 -0.0032200 -1.9e-06 0.0001720 NA NA NA 10 2222.169 -4424.249 38.93904 3.497280e-09
9 -2.292576 -7.9e-06 0.0006068 -1e-07 -0.0006129 0.0004487 0.0012922 -1.0e-06 NA NA 0.0003117 NA 10 2222.151 -4424.212 38.97548 3.434136e-09
20 -2.292081 -7.9e-06 0.0005863 -1e-07 -0.0006079 0.0003718 0.0020931 NA NA -1e-07 NA -0.0002167 10 2219.415 -4418.741 44.44700 2.226858e-10
26 -2.292680 -7.7e-06 0.0006054 -1e-07 -0.0006090 0.0003997 -0.0016523 NA 0.0000905 NA 0.0003035 -0.0002108 11 2220.273 -4418.437 44.75033 1.913486e-10
12 -2.292870 -8.0e-06 0.0006222 -1e-07 -0.0006246 0.0004485 0.0023252 NA NA -1e-07 0.0003115 NA 10 2219.171 -4418.253 44.93462 1.745048e-10
6 -2.291537 -7.8e-06 0.0005690 -1e-07 -0.0005789 0.0003970 -0.0002205 NA 0.0000895 -1e-07 NA NA 10 2219.115 -4418.140 45.04770 1.649120e-10
19 -2.291677 -7.8e-06 0.0005670 -1e-07 -0.0005823 0.0003628 -0.0022333 -1.1e-06 0.0001231 NA NA -0.0001991 11 2215.992 -4409.877 53.31107 2.647794e-12
25 -2.292858 -7.8e-06 0.0006128 -1e-07 -0.0006237 0.0004086 0.0009470 -3.0e-07 NA NA 0.0003030 -0.0002218 11 2215.932 -4409.756 53.43142 2.493164e-12
5 -2.292014 -8.2e-06 0.0005902 -1e-07 -0.0006074 0.0004143 0.0028381 -1.4e-06 NA -2e-07 NA NA 10 2214.737 -4409.384 53.80399 2.069413e-12
11 -2.292555 -8.0e-06 0.0006063 -1e-07 -0.0006006 0.0004350 -0.0029759 -1.7e-06 0.0001609 NA 0.0003077 NA 11 2215.650 -4409.191 53.99633 1.879667e-12
28 -2.293197 -8.0e-06 0.0006311 -1e-07 -0.0006382 0.0004180 0.0021203 NA NA -1e-07 0.0003025 -0.0001829 11 2212.889 -4403.671 59.51672 1.189444e-13
22 -2.291971 -7.9e-06 0.0005818 -1e-07 -0.0005983 0.0003671 0.0002469 NA 0.0000652 -1e-07 NA -0.0002034 11 2212.849 -4403.591 59.59667 1.142835e-13
14 -2.292756 -7.9e-06 0.0006172 -1e-07 -0.0006133 0.0004393 -0.0001273 NA 0.0000859 -1e-07 0.0003107 NA 11 2212.598 -4403.089 60.09851 8.892208e-14
21 -2.292218 -8.1e-06 0.0005934 -1e-07 -0.0006146 0.0003809 0.0024002 -8.0e-07 NA -1e-07 NA -0.0001865 11 2208.522 -4394.936 68.25187 1.508452e-15
27 -2.292778 -7.9e-06 0.0006110 -1e-07 -0.0006121 0.0004077 -0.0021330 -9.0e-07 0.0001191 NA 0.0003020 -0.0001710 12 2209.466 -4394.805 68.38296 1.412746e-15
7 -2.291971 -8.2e-06 0.0005884 -1e-07 -0.0005954 0.0004025 -0.0009542 -2.0e-06 0.0001395 -1e-07 NA NA 11 2208.232 -4394.357 68.83056 1.129455e-15
13 -2.293157 -8.2e-06 0.0006351 -1e-07 -0.0006382 0.0004545 0.0027477 -1.2e-06 NA -2e-07 0.0003072 NA 11 2208.216 -4394.325 68.86300 1.111284e-15
30 -2.293086 -8.0e-06 0.0006266 -1e-07 -0.0006285 0.0004133 0.0002597 NA 0.0000657 -1e-07 0.0003026 -0.0001694 12 2206.324 -4388.520 74.66758 6.100666e-17
23 -2.292138 -8.1e-06 0.0005912 -1e-07 -0.0006038 0.0003798 -0.0003934 -1.4e-06 0.0001065 -1e-07 NA -0.0001426 12 2202.054 -4379.981 83.20700 8.532288e-19
29 -2.293309 -8.1e-06 0.0006369 -1e-07 -0.0006437 0.0004255 0.0023796 -7.0e-07 NA -1e-07 0.0003014 -0.0001575 12 2201.995 -4379.862 83.32553 8.041324e-19
15 -2.293107 -8.3e-06 0.0006330 -1e-07 -0.0006269 0.0004432 -0.0007647 -1.7e-06 0.0001292 -1e-07 0.0003042 NA 12 2201.708 -4379.289 83.89845 6.038371e-19
31 -2.293228 -8.2e-06 0.0006348 -1e-07 -0.0006332 0.0004243 -0.0003138 -1.2e-06 0.0001027 -1e-07 0.0003006 -0.0001152 13 2195.526 -4364.904 98.28412 4.540571e-22
#Run model for DR
Subset.MCC.DR.top <- gls(MCC.DR ~ Sexual_selection_ppca
                         + log(range.size.m2)
                         + bioclim4 #Seasonal variation
                         + residuals.PC1 #Spatial variation
                         + PC1.LIG #Long-term climate variation
                         + NPP,
                correlation = corPagel(1, phy = pruned.MCC.Subset.tree, fixed = FALSE), 
                data = SS.subset, 
                method = "REML")
saveRDS(Subset.MCC.DR.top, 'data/Subset.MCC.DR.top.rds')

#Run model for ND
Subset.MCC.ND.top <- gls(MCC.ND ~ Sexual_selection_ppca
                         + log(range.size.m2)
                         + bioclim4 #Seasonal variation
                         + residuals.PC1 #Spatial variation
                         + PC1.LIG #Long-term climate variation
                         + NPP,
                correlation = corPagel(1, phy = pruned.MCC.Subset.tree, fixed = TRUE), #lambda = 1
                data = SS.subset, 
                method = "REML")
saveRDS(Subset.MCC.ND.top, 'data/Subset.MCC.ND.top.rds')

#Run the 100 models for DR and ND using the best model:
Subset.data.noMCC <- SS.subset %>% dplyr::select(TipLabel, Sexual_selection_ppca, range.size.m2, bioclim4, residuals.PC1, PC1.LIG, NPP)

#Take the restricted data and make it simpler with just responses and predictors.Note that we join the es.values for the 100 trees
Subset.DR.model.data <- lapply(es.list, function(x) { #es.list is a list of ES values calculated earlier
  left_join(Subset.data.noMCC, 
            x %>% as.data.frame() %>% tibble::rownames_to_column() %>% `colnames<-`(c("TipLabel", "DR")), 
            by = "TipLabel")
})

#PGLS needs tiplabel as rowname
Subset.DR.model.data <- lapply(Subset.DR.model.data, function(x) {
  tibble::column_to_rownames(x, "TipLabel")})

#Prune the trees
Subset.pruned.trees<-lapply(passerine.trees, function(x) {
  drop.tip(x,x$tip.label[-match(SS.subset$TipLabel, x$tip.label)])
})

#Use mapply to create a list of PGLS global models
Subset.DR.pgls.models <- mcmapply(function(x,y) {
  gls(DR ~ Sexual_selection_ppca 
         + log(range.size.m2)
         + bioclim4 #Seasonal variation
         + residuals.PC1 #Spatial variation
         + PC1.LIG #Long-term climate variation
         + NPP,
    correlation = corPagel(0.9711, phy = y, fixed = TRUE), 
    data = x, 
    method = "REML")
}, x = Subset.DR.model.data, y = Subset.pruned.trees,
SIMPLIFY = FALSE, #Prevents the catastrophic loss of structure, names and attributes
mc.cores = 8) #Specify core number 

saveRDS(Subset.DR.pgls.models, "data/Subset.DR.pgls.models.rds")

#Now for Node Density:
Subset.ND.model.data <- lapply(nd.list, function(x) { #es.list is a list of ES values calculated earlier
  left_join(Subset.data.noMCC, 
            x %>% as.data.frame() %>% tibble::rownames_to_column() %>% `colnames<-`(c("TipLabel", "ND")), 
            by = "TipLabel")
})

#PGLS needs tiplabel as rowname
Subset.ND.model.data <- lapply(Subset.ND.model.data, function(x) {
  tibble::column_to_rownames(x, "TipLabel")})

#Use mapply to create a list of PGLS global models
Subset.ND.pgls.models <- mcmapply(function(x,y) {
gls(ND ~ Sexual_selection_ppca 
         + log(range.size.m2)
         + bioclim4 #Seasonal variation
         + residuals.PC1 #Spatial variation
         + PC1.LIG #Long-term climate variation
         + NPP,
    corPagel(1, phy = y, fixed = TRUE), 
    data = x, 
    method = "REML")
}, x = Subset.ND.model.data, y = Subset.pruned.trees, 
SIMPLIFY = FALSE, #Prevents the catastrophic loss of structure, names and attributes
mc.cores = 8) #Specify core number 

saveRDS(Subset.ND.pgls.models, "data/Subset.ND.pgls.models.rds")

#BAMM Top Models 

Subset.MCC.Lambda.top <- gls(mean.lambda ~ Sexual_selection_ppca
                      + log(range.size.m2)
                      + bioclim4 #Seasonal variation
                      + residuals.PC1 #Spatial variation
                      + PC1.LIG #Long-term climate variation
                      + NPP,
                weights = ~ sqrt(var.lambda),
                correlation = corBrownian(phy = pruned.MCC.Subset.tree), 
                data = SS.subset, 
                method = "REML")
saveRDS(Subset.MCC.Lambda.top, 'data/Subset.MCC.Lambda.top.rds')

Subset.MCC.Mu.top <- gls(mean.mu ~ Sexual_selection_ppca
                      + log(range.size.m2)
                      + bioclim4 #Seasonal variation
                      + residuals.PC1 #Spatial variation
                      + PC1.LIG #Long-term climate variation
                      + NPP,
                weights = ~ sqrt(var.mu),
                correlation = corBrownian(phy = pruned.MCC.Subset.tree), 
                data = SS.subset, 
                method = "REML")
saveRDS(Subset.MCC.Mu.top, 'data/Subset.MCC.Mu.top.rds')

#Now for BAMM 100 models

Subset.BAMM.model.data <- lapply(BAMM.df, function(x) { #es.list is a list of ES values calculated earlier
  left_join(Subset.data.noMCC %>% dplyr::select(TipLabel, Sexual_selection_ppca, range.size.m2, bioclim4, residuals.PC1, PC1.LIG, NPP), 
            x %>% as.data.frame(), 
            by = "TipLabel")
})

#PGLS needs tiplabel as rowname
Subset.BAMM.model.data <- lapply(Subset.BAMM.model.data, function(x) {
  tibble::column_to_rownames(x, "TipLabel")})

#Use mapply to create a list of PGLS global models
Subset.BAMM.lambda.pgls.models <- mcmapply(function(x,y) {
  gls(mean.lambda ~ Sexual_selection_ppca 
         + log(range.size.m2)
         + bioclim4 #Seasonal variation
         + residuals.PC1 #Spatial variation
         + PC1.LIG #Long-term climate variation
         + NPP,
    weights = ~ sqrt(var.lambda),      
    corBrownian(phy = y), #lambda = 1
    data = x, 
    method = "REML")
}, x = Subset.BAMM.model.data, y = Subset.pruned.trees,
SIMPLIFY = FALSE, #Prevents the catastrophic loss of structure, names and attributes
mc.cores = 8) #Specify core number 

saveRDS(Subset.BAMM.lambda.pgls.models, "data/Subset.BAMM.lambda.pgls.models.rds")

#Use mapply to create a list of PGLS global models
Subset.BAMM.mu.pgls.models <- mcmapply(function(x,y) {
  gls(mean.mu ~ Sexual_selection_ppca
         + log(range.size.m2)
         + bioclim4 #Seasonal variation
         + residuals.PC1 #Spatial variation
         + PC1.LIG #Long-term climate variation
         + NPP,
    weights = ~ sqrt(var.mu),      
    corBrownian(phy = y), #lambda = 1
    data = x, 
    method = "REML")
}, x = Subset.BAMM.model.data, y = Subset.pruned.trees,
SIMPLIFY = FALSE, #Prevents the catastrophic loss of structure, names and attributes
mc.cores = 8) #Specify core number 

saveRDS(Subset.BAMM.mu.pgls.models, "data/Subset.BAMM.mu.pgls.models.rds")

The phylogenetic signal from the DR model is ~0.97. The other models had phylogenetic signals ~ 1, as such they were fixed at 1 to avoid problems of convergence.

Subset.MCC.DR.top <- readRDS('data/Subset.MCC.DR.top.rds')
Subset.MCC.DR.top[["modelStruct"]][["corStruct"]] %>% `names<-`("DR lambda") %>% pander()
  • DR lambda: 0.9711
Subset.DR.pgls.models <- readRDS("data/Subset.DR.pgls.models.rds")
Subset.MCC.DR.top <- readRDS('data/Subset.MCC.DR.top.rds')

Subset.cols <- brewer.pal(n = 7, name = "Dark2")[-6]

Subset.DR.pgls.summary <- lapply(Subset.DR.pgls.models, function(x) {
  data.frame(x$coefficients, confint(x)) %>% tibble::rownames_to_column()
})

Subset.MCC.DR.summary <- data.frame(Subset.MCC.DR.top$coefficients, confint(Subset.MCC.DR.top)) %>% tibble::rownames_to_column()

Subset.DR.pgls.summary <- bind_rows(Subset.DR.pgls.summary) %>% `colnames<-`(c("Parameter", "Estimate", "LCI", "UCI"))

colnames(Subset.MCC.DR.summary) <- c("Parameter", "Estimate", "LCI", "UCI")

parameter_names <- c(
                    `bioclim4` = "Temperature Seasonality",
                    `log(range.size.m2)` = "Range Size (log-transformed)",
                    `NPP` = "NPP",
                    `PC1.LIG` = "Long-term Temperature Variation",
                    `residuals.PC1` = "Spatial Temperature Variation",
                    `Sexual_selection_ppca` = "Sexual Selection"
                    )

Subset.DR.pgls.summary$Parameter = factor(Subset.DR.pgls.summary$Parameter, levels=c('bioclim4','log(range.size.m2)','NPP','PC1.LIG', 'residuals.PC1', 'Sexual_selection_ppca'))

Subset.MCC.DR.summary$Parameter = factor(Subset.MCC.DR.summary$Parameter, levels=c('bioclim4','log(range.size.m2)','NPP','PC1.LIG', 'residuals.PC1', 'Sexual_selection_ppca'))

Subset.DR.plot <-Subset.DR.pgls.summary %>% filter(Parameter != "(Intercept)") %>% 
  ggplot(aes(x = Estimate, y = Parameter, fill = Parameter)) +
  geom_density(aes(y = ..scaled..),
               position = position_nudge(y= 0.6))+
  geom_jitter(shape = 21, alpha = 0.5, size = 0.75)+
  geom_vline(xintercept = 0)+
  theme_minimal()+
  theme(axis.text.y = element_blank(),
        legend.position = "none")

Subset.DR.plot <- Subset.DR.plot + geom_errorbarh(data = Subset.MCC.DR.summary %>% filter(Parameter != "(Intercept)"), 
                 aes(xmin = LCI,
                     xmax = UCI, y = Parameter,
                     color = Parameter), 
                 height = 0, show.legend = F, 
                 position = position_nudge(y= -0.75))+
  
  geom_point(data = Subset.MCC.DR.summary %>% filter(Parameter != "(Intercept)"), 
             aes(x = Estimate, y = Parameter, fill = Parameter), 
             shape=21, color = "grey20",
             size = 3,
             position = position_nudge(y= -0.75)) +
  scale_y_discrete(expand = c(0.5,.5))+
  facet_wrap(~ Parameter, scales = "free", nrow = 6, labeller = as_labeller(parameter_names))+
  scale_fill_manual(values = Subset.cols)+
  scale_color_manual(values = Subset.cols)+
  theme(panel.grid.major.y = element_blank(),
        axis.title.y = element_blank(),
        axis.text.x = element_text(size = 7),
        plot.margin=unit(c(.5,.5,.5,.5),"cm"))+
  ggtitle("Male-Bias SS DR")

####################################
#For ND

Subset.ND.pgls.models <- readRDS("data/Subset.ND.pgls.models.rds")
Subset.MCC.ND.top <- readRDS('data/Subset.MCC.ND.top.rds')

Subset.ND.pgls.summary <- lapply(Subset.ND.pgls.models, function(x) {
  data.frame(x$coefficients, confint(x)) %>% tibble::rownames_to_column()
})

Subset.MCC.ND.summary <- data.frame(Subset.MCC.ND.top$coefficients, confint(Subset.MCC.ND.top)) %>% tibble::rownames_to_column()

Subset.ND.pgls.summary <- bind_rows(Subset.ND.pgls.summary) %>% `colnames<-`(c("Parameter", "Estimate", "LCI", "UCI"))

colnames(Subset.MCC.ND.summary) <- c("Parameter", "Estimate", "LCI", "UCI")

Subset.ND.pgls.summary$Parameter = factor(Subset.ND.pgls.summary$Parameter, levels=c('bioclim4','log(range.size.m2)','NPP','PC1.LIG', 'residuals.PC1', 'Sexual_selection_ppca'))

Subset.MCC.ND.summary$Parameter = factor(Subset.MCC.ND.summary$Parameter, levels=c('bioclim4','log(range.size.m2)','NPP','PC1.LIG', 'residuals.PC1', 'Sexual_selection_ppca'))

Subset.ND.plot <-Subset.ND.pgls.summary %>% filter(Parameter != "(Intercept)") %>% 
  ggplot(aes(x = Estimate, y = Parameter, fill = Parameter)) +
  geom_density(aes(y = ..scaled..),
               position = position_nudge(y= 0.6))+
  geom_jitter(shape = 21, alpha = 0.5, size = 0.75)+
  geom_vline(xintercept = 0)+
  theme_minimal()+
  theme(axis.text.y = element_blank(),
        legend.position = "none")

Subset.ND.plot <- Subset.ND.plot + geom_errorbarh(data = Subset.MCC.ND.summary %>% filter(Parameter != "(Intercept)"), 
                 aes(xmin = LCI,
                     xmax = UCI, y = Parameter,
                     color = Parameter), 
                 height = 0, show.legend = F, 
                 position = position_nudge(y= -0.75))+
  
  geom_point(data = Subset.MCC.ND.summary %>% filter(Parameter != "(Intercept)"), 
             aes(x = Estimate, y = Parameter, fill = Parameter), 
             shape=21, color = "grey20",
             size = 3,
             position = position_nudge(y= -0.75)) +
  scale_y_discrete(expand = c(0.5,.5))+
  facet_wrap(~ Parameter, scales = "free", nrow = 6, labeller = as_labeller(parameter_names))+
  scale_fill_manual(values = Subset.cols)+
  scale_color_manual(values = Subset.cols)+
  theme(panel.grid.major.y = element_blank(),
        axis.title.y = element_blank(),
        axis.text.x = element_text(size = 7),
        plot.margin=unit(c(.5,.5,.5,.5),"cm"))+
  ggtitle("Male-Bias SS ND")

######################

#For Lambda
Subset.BAMM.lambda.pgls.models <- readRDS("data/Subset.BAMM.lambda.pgls.models.rds")
Subset.MCC.Lambda.top <- readRDS('data/Subset.MCC.Lambda.top.rds')

Subset.Lambda.pgls.summary <- lapply(Subset.BAMM.lambda.pgls.models, function(x) {
  data.frame(x$coefficients, confint(x)) %>% tibble::rownames_to_column()
})

Subset.MCC.Lambda.summary <- data.frame(Subset.MCC.Lambda.top$coefficients, confint(Subset.MCC.Lambda.top)) %>% tibble::rownames_to_column()

Subset.Lambda.pgls.summary <- bind_rows(Subset.Lambda.pgls.summary) %>% `colnames<-`(c("Parameter", "Estimate", "LCI", "UCI"))

colnames(Subset.MCC.Lambda.summary) <- c("Parameter", "Estimate", "LCI", "UCI")

Subset.Lambda.pgls.summary$Parameter = factor(Subset.Lambda.pgls.summary$Parameter, levels=c('bioclim4','log(range.size.m2)','NPP','PC1.LIG', 'residuals.PC1', 'Sexual_selection_ppca'))

Subset.MCC.Lambda.summary$Parameter = factor(Subset.MCC.Lambda.summary$Parameter, levels=c('bioclim4','log(range.size.m2)','NPP','PC1.LIG', 'residuals.PC1', 'Sexual_selection_ppca'))

Subset.Lambda.pgls.summary.RO <- dcast(Subset.Lambda.pgls.summary %>% filter(Parameter != "(Intercept)"), Estimate ~ Parameter, value.var = "Estimate")
Subset.Lambda.pgls.summary.RO$Estimate <- NULL
Subset.Lambda.pgls.summary.RO <- sapply(Subset.Lambda.pgls.summary.RO, function(x) {
  remove_outliers(x, na.rm = T)})
Subset.Lambda.pgls.summary.RO <-melt(Subset.Lambda.pgls.summary.RO) %>% na.omit()
Subset.Lambda.pgls.summary.RO$Var1 <- NULL
colnames(Subset.Lambda.pgls.summary.RO) <- c("Parameter", "Estimate")


Subset.Lambda.plot <-Subset.Lambda.pgls.summary.RO %>% 
  ggplot(aes(x = Estimate, y = Parameter, fill = Parameter)) +
  geom_density(aes(y = ..scaled..),
               position = position_nudge(y= 0.6))+
  geom_jitter(shape = 21, alpha = 0.5, size = 0.75)+
  geom_vline(xintercept = 0)+
  theme_minimal()+
  theme(axis.text.y = element_blank(),
        legend.position = "none")

Subset.Lambda.plot <- Subset.Lambda.plot + geom_errorbarh(data = Subset.MCC.Lambda.summary %>% filter(Parameter != "(Intercept)"), 
                 aes(xmin = LCI,
                     xmax = UCI, y = Parameter,
                     color = Parameter), 
                 height = 0, show.legend = F, 
                 position = position_nudge(y= -0.75))+
  
  geom_point(data = Subset.MCC.Lambda.summary %>% filter(Parameter != "(Intercept)"), 
             aes(x = Estimate, y = Parameter, fill = Parameter), 
             shape=21, color = "grey20",
             size = 3,
             position = position_nudge(y= -0.75)) +
  scale_y_discrete(expand = c(0.5,.5))+
  facet_wrap(~ Parameter, scales = "free", nrow = 6, labeller = as_labeller(parameter_names))+
  scale_fill_manual(values = Subset.cols)+
  scale_color_manual(values = Subset.cols)+
  theme(panel.grid.major.y = element_blank(),
        axis.title.y = element_blank(),
        axis.text.x = element_text(size = 7),
        plot.margin=unit(c(.5,.5,.5,.5),"cm"))+
  ggtitle("Male-Bias SS Speciaton")

######################

#For Mu

Subset.BAMM.Mu.pgls.models <- readRDS("data/Subset.BAMM.mu.pgls.models.rds")
Subset.MCC.Mu.top <- readRDS('data/Subset.MCC.Mu.top.rds')

Subset.Mu.pgls.summary <- lapply(Subset.BAMM.Mu.pgls.models, function(x) {
  data.frame(x$coefficients, confint(x)) %>% tibble::rownames_to_column()
})

Subset.MCC.Mu.summary <- data.frame(Subset.MCC.Mu.top$coefficients, confint(Subset.MCC.Mu.top)) %>% tibble::rownames_to_column()

Subset.Mu.pgls.summary <- bind_rows(Subset.Mu.pgls.summary) %>% `colnames<-`(c("Parameter", "Estimate", "LCI", "UCI"))

colnames(Subset.MCC.Mu.summary) <- c("Parameter", "Estimate", "LCI", "UCI")

Subset.Mu.pgls.summary$Parameter = factor(Subset.Mu.pgls.summary$Parameter, levels=c('bioclim4','log(range.size.m2)','NPP','PC1.LIG', 'residuals.PC1', 'Sexual_selection_ppca'))

Subset.MCC.Mu.summary$Parameter = factor(Subset.MCC.Mu.summary$Parameter, levels=c('bioclim4','log(range.size.m2)','NPP','PC1.LIG', 'residuals.PC1', 'Sexual_selection_ppca'))


Subset.Mu.pgls.summary.RO <- dcast(Subset.Mu.pgls.summary %>% filter(Parameter != "(Intercept)"), Estimate ~ Parameter, value.var = "Estimate")
Subset.Mu.pgls.summary.RO$Estimate <- NULL
Subset.Mu.pgls.summary.RO <- sapply(Subset.Mu.pgls.summary.RO, function(x) {
  remove_outliers(x, na.rm = T)})
Subset.Mu.pgls.summary.RO <-melt(Subset.Mu.pgls.summary.RO) %>% na.omit()
Subset.Mu.pgls.summary.RO$Var1 <- NULL
colnames(Subset.Mu.pgls.summary.RO) <- c("Parameter", "Estimate")

Subset.Mu.plot <-Subset.Mu.pgls.summary.RO %>% 
  ggplot(aes(x = Estimate, y = Parameter, fill = Parameter)) +
  geom_density(aes(y = ..scaled..),
               position = position_nudge(y= 0.6))+
  geom_jitter(shape = 21, alpha = 0.5, size = 0.75)+
  geom_vline(xintercept = 0)+
  theme_minimal()+
  theme(axis.text.y = element_blank(),
        legend.position = "none")

Subset.Mu.plot <- Subset.Mu.plot + geom_errorbarh(data = Subset.MCC.Mu.summary %>% filter(Parameter != "(Intercept)"), 
                 aes(xmin = LCI,
                     xmax = UCI, y = Parameter,
                     color = Parameter), 
                 height = 0, show.legend = F, 
                 position = position_nudge(y= -0.75))+
  
  geom_point(data = Subset.MCC.Mu.summary %>% filter(Parameter != "(Intercept)"), 
             aes(x = Estimate, y = Parameter, fill = Parameter), 
             shape=21, color = "grey20",
             size = 3,
             position = position_nudge(y= -0.75)) +
  scale_y_discrete(expand = c(0.5,.5))+
  facet_wrap(~ Parameter, scales = "free", nrow = 6, labeller = as_labeller(parameter_names))+
  scale_fill_manual(values = Subset.cols)+
  scale_color_manual(values = Subset.cols)+
  theme(panel.grid.major.y = element_blank(),
        axis.title.y = element_blank(),
        axis.text.x = element_text(size = 7),
        plot.margin=unit(c(.5,.5,.5,.5),"cm"))+
  ggtitle("Male-Bias SS Extinction")

#########################################
grid.arrange(symmetrise_scale(Subset.DR.plot, "x"),
             symmetrise_scale(Subset.ND.plot, "x"),
             symmetrise_scale(Subset.Lambda.plot, "x"), 
             symmetrise_scale(Subset.Mu.plot, "x"), 
             nrow = 1)

Figure S12: Measures of sexual selection a phylogenetic pca of sexual dimorphism, social polygyny and paternal care provide some evidence to suggest sexual selection predicts extinction rate (DR and ND models). This figure is seen in Figure 1 within the manuscript

Table S14: MCC model estimates from the above model are presented here as numberical values with 95 % CIs.

Subset.MCC.DR.summary %>% filter(Parameter != "(Intercept)") %>% pander(caption = "MCC DR Estimates from Male-bias SS")
MCC DR Estimates from Male-bias SS
Parameter Estimate LCI UCI
Sexual_selection_ppca 0.03887 0.009472 0.06827
log(range.size.m2) -0.01131 -0.02055 -0.002058
bioclim4 2.448e-05 -6.115e-05 0.0001101
residuals.PC1 0.0004409 -0.00859 0.009472
PC1.LIG 0.006881 -0.001666 0.01543
NPP -7.335e-07 -3.45e-06 1.984e-06
Subset.MCC.ND.summary %>% filter(Parameter != "(Intercept)") %>% pander(caption = "MCC ND Estimates from Male-bias SS")
MCC ND Estimates from Male-bias SS
Parameter Estimate LCI UCI
Sexual_selection_ppca 0.0004381 -0.0004178 0.001294
log(range.size.m2) -0.00018 -0.0004838 0.0001238
bioclim4 -4.586e-07 -3.601e-06 2.683e-06
residuals.PC1 7.875e-05 -0.0002533 0.0004108
PC1.LIG 0.0001073 -0.0001577 0.0003723
NPP -1.723e-08 -1.229e-07 8.84e-08
Subset.MCC.Lambda.summary %>% filter(Parameter != "(Intercept)") %>% pander(caption = "MCC BAMM Speciation Estimates from Male-bias SS")
MCC BAMM Speciation Estimates from Male-bias SS
Parameter Estimate LCI UCI
Sexual_selection_ppca 0.0009422 -0.002845 0.004729
log(range.size.m2) 0.0005471 -0.0006171 0.001711
bioclim4 -7.604e-06 -1.939e-05 4.177e-06
residuals.PC1 0.0004011 -0.0008461 0.001648
PC1.LIG -0.000567 -0.001625 0.0004906
NPP -5.98e-08 -4.751e-07 3.555e-07
Subset.MCC.Mu.summary %>% filter(Parameter != "(Intercept)") %>% pander(caption = "MCC BAMM Extinction Estimates from Male-bias SS")
MCC BAMM Extinction Estimates from Male-bias SS
Parameter Estimate LCI UCI
Sexual_selection_ppca 0.003443 -0.004091 0.01098
log(range.size.m2) 0.0004388 -0.002251 0.003129
bioclim4 -1.361e-05 -4.164e-05 1.442e-05
residuals.PC1 0.000148 -0.002828 0.003124
PC1.LIG -0.0009247 -0.003254 0.001405
NPP -2.112e-07 -1.169e-06 7.471e-07

Table S15: Highest posterior density intervals are calculated for the above figure (ie; models using Male-bias sexual selection measures). using this measure the models using \(\lambda_{DR}\) have model estimates that fall on the positive side 95 % of the time in 100 trees. For \(\lambda_{ND}\), there is a positive skew, however the 95 % HPD interval overlaps zero. \(\lambda_{BAMM}\) also shows a lesser positive skew. The HPD range is determined using the hdi function with a 95 % credible interval. These intervals do not take into account the variance associated with each interval and thus are not an estimate of model precision. Intervals not overlapping zero suggest that 95 % of trees from the posterior generate a model estimate for the given parameter that in the same direction (+ or -). These intervals are calculated in the same way as in Table S7 and Table S9.

Subset.hpd.DR.top <- list()
Subset.DR.pgls.summary <- na.omit(Subset.DR.pgls.summary)
for(x in unique(Subset.DR.pgls.summary$Parameter)) {
Subset.hpd.DR.top[[x]] = hdi(Subset.DR.pgls.summary %>% filter(Parameter == x) %>% dplyr::select(Estimate))
}
Subset.hpd.DR.top <- bind_rows(Subset.hpd.DR.top) %>% `rownames<-`(c("Lower", "Upper"))


saveRDS(Subset.hpd.DR.top, 'data/Subset.hpd.DR.top.rds')

#For ND
Subset.hpd.ND.top <- list()
Subset.ND.pgls.summary <- na.omit(Subset.ND.pgls.summary)
for(x in unique(Subset.ND.pgls.summary$Parameter)) {
Subset.hpd.ND.top[[x]] = hdi(Subset.ND.pgls.summary %>% filter(Parameter == x) %>% dplyr::select(Estimate))
}
Subset.hpd.ND.top <- bind_rows(Subset.hpd.ND.top) %>% `rownames<-`(c("Lower", "Upper"))

saveRDS(Subset.hpd.ND.top, 'data/Subset.hpd.ND.top.rds')

###############
Subset.hpd.Lambda.top <- list()
Subset.Lambda.pgls.summary <- na.omit(Subset.Lambda.pgls.summary)
for(x in unique(Subset.Lambda.pgls.summary$Parameter)) {
Subset.hpd.Lambda.top[[x]] = hdi(Subset.Lambda.pgls.summary %>% filter(Parameter == x) %>% dplyr::select(Estimate))
}
Subset.hpd.Lambda.top <- bind_rows(Subset.hpd.Lambda.top) %>% `rownames<-`(c("Lower", "Upper")) 
saveRDS(Subset.hpd.Lambda.top, 'data/Subset.hpd.Lambda.top.rds')

################
Subset.hpd.Mu.top <- list()
Subset.Mu.pgls.summary <- na.omit(Subset.Mu.pgls.summary)
for(x in unique(Subset.Mu.pgls.summary$Parameter)) {
Subset.hpd.Mu.top[[x]] = hdi(Subset.Mu.pgls.summary %>% filter(Parameter == x) %>% dplyr::select(Estimate))
}
Subset.hpd.Mu.top <- bind_rows(Subset.hpd.Mu.top) %>% `rownames<-`(c("Lower", "Upper"))

saveRDS(Subset.hpd.Mu.top, 'data/Subset.hpd.Mu.top.rds')

Subset.hpd.DR.top %>% pander(split.table = Inf, digits = 3, caption = "Subset DR HPD Interval")
Subset DR HPD Interval
  Sexual_selection_ppca log(range.size.m2) bioclim4 residuals.PC1 PC1.LIG NPP
Lower 0.00911 -0.0158 -9.31e-05 0.000313 -0.0036 -2.96e-06
Upper 0.0608 -0.000501 2.76e-05 0.0128 0.00894 1.19e-06
Subset.hpd.ND.top %>% pander(split.table = Inf, digits = 3, caption = "Subset ND HPD Interval")
Subset ND HPD Interval
  Sexual_selection_ppca log(range.size.m2) bioclim4 residuals.PC1 PC1.LIG NPP
Lower -0.00031 -0.000426 -3.63e-06 -3.05e-05 -0.000249 -9.23e-08
Upper 0.00154 0.000166 1.64e-06 0.000487 0.000249 6.61e-08
Subset.hpd.Lambda.top %>% pander(split.table = Inf, digits = 3, caption = "Subset BAMM Speciation HPD Interval")
Subset BAMM Speciation HPD Interval
  Sexual_selection_ppca log(range.size.m2) bioclim4 residuals.PC1 PC1.LIG NPP
Lower -0.013 -0.0126 -9.69e-05 -0.00745 -0.00989 -2.09e-06
Upper 0.0309 0.0166 5.42e-05 0.0187 0.0047 1.99e-06
Subset.hpd.Mu.top %>% pander(split.table = Inf, digits = 3, caption = "Subset BAMM Extinction HPD Interval")
Subset BAMM Extinction HPD Interval
  Sexual_selection_ppca log(range.size.m2) bioclim4 residuals.PC1 PC1.LIG NPP
Lower -0.0516 -0.0222 -0.000423 -0.0328 -0.0186 -3.14e-06
Upper 0.0576 0.0306 0.000171 0.0207 0.0158 2.78e-06

Phylogenetic Path Analysis

We undertook a phylogenetic path analysis using a model of causal pathways a priori. To do this we use the phylopath function. Our path analysis is undertaken on the dataset subset to species with measures of male bias sexual selection. Additionallty, the path analysis is restricted to using one measure of speciation (\(\lambda_{DR}\)) and the MCC tree. Our reasoning for the causal links are described breifly below:

  1. \(\lambda_{DR}\) dependent upon Sexual Dichromatism: Sexual dichromatism would directly impact speciation if there is an association with the evolvability of chromatism and speciation, this may more rapidly lead to niche divergence or reproductive isolation.
  2. \(\lambda_{DR}\) dependent upon Sexual Selection: This was a major hypothesis in our study (see introduction).
  3. \(\lambda_{DR}\) dependent upon Temperature Seasonality: More variable environments are broadly expected to impact the resources, niches and fitness of species; thus there is an expectation that more variable environments may impact speciation rate.
  4. \(\lambda_{DR}\) dependent upon Range Size: As previously seen in the PGLS models range size shows a negative correlation with speciation rate; expected reasons for this causal link include lower dispersal ability in smaller ranges correlating with lower gene flow and greater speciation. It may also reflect niche specialisation and thus in line with the previous point increase speciation and reduce gene flow.
  5. Sexual Dichromatism dependent upon Sexual Selection: Higher levels of sexual selection will promote showy males and drab females (sexual dichromatism).
  6. Sexual Dichromatism dependent upon Temperature Seasonality: Ecological aspects may facilitate sexual dimorphism as niche partitioning between the sexes causes sex-limited traits and colouration that are independent of mate choice or male-male competition.
  7. Sexual Selection dependent upon Temperature Seasonality: Greater environmental variability/stress may impact the usefulness (fitness effects) of sexual selection (see introduction).
  8. Range Size dependent upon Sexual Selection: Sexual selection may often depend on local adaptations and the alignment of sexually selected traits with locally valuable traits, range size may limit this leading to a negative correlation. Alternatively if sexual selection has added fitness and evolvability benefits there may be increased range sizes in sexually selected species.
  9. Range Size dependent upon Temperature Seasonality: Anticipated to be more of a correlation than causal link. Yet a speecies range size may be dependent upon how variable their environment is and the availability to local resources.
#Set rownames to match tree
rownames(SS.subset) <- SS.subset$TipLabel
pruned.MCC.Subset.tree <- readRDS('data/pruned.MCC.Subset.tree.rds')
SS.subset$log.range.size <- log(SS.subset$range.size.m2)

SS.subset2 <- SS.subset %>% rename(
  DR = MCC.DR,
  SD = SDi,
  TS = bioclim4,
  RS = log.range.size,
  SS = Sexual_selection_ppca
)

models.Subset <- define_model_set(
  one = c(DR ~ SD, 
          DR ~ SS,
          DR ~ TS,
          DR ~ RS,
          SD ~ SS,
          SD ~ TS,
          SS ~ TS,
          RS ~ SS,
          RS ~ TS)
)

result <- phylo_path(models.Subset, data = SS.subset2, tree = pruned.MCC.Subset.tree, model = 'lambda')
  
#best_model <- best(result, boot = 500)
set.seed(1)
path.plot <- plot(x = best_model,
     type = "color",
     algorithm = 'gem',
     manual_layout = NULL,
     curvature = 0.1,
     colors = c("#b2182b", "#2166ac"),
     show.legend = F)

saveRDS(best_model, 'data/path_model.rds')
#Inspect the coefficients with their SE

# coef_plot(best_model, error_bar = "ci", reverse_order = TRUE) + 
#   ggplot2::coord_flip()+
#   ggplot2::theme_minimal()

# pdf("Figures/Path_Plot.pdf", width=8, height=8)
# path.plot
# dev.off()

# path.plot
best_model <- readRDS('data/path_model.rds')
coef_plot(best_model, error_bar = "ci", reverse_order = TRUE) + 
  ggplot2::coord_flip()+
  ggplot2::theme_minimal()

Figure S13: Path analysis standardized regression coefficients vary across relationships. Error bars are derived from confidence intervals through 500 bootstrapped iterations. Paths with error bars not everlapping zero are presented with an asterix in Figure 3 within the manuscript.

Additional Figures and Tables

BAMM.df <- readRDS('data/BAMM.df.rds')
BAMM.Rates.100Trees <- do.call(rbind, BAMM.df)
BAMM.Rates.100Trees$mean.lambda <- BAMM.Rates.100Trees$mean.lambda
BAMM.Lambda.100Trees <- ggiraphExtra::summarySE(BAMM.Rates.100Trees, measurevar = "mean.lambda", groupvars = "TipLabel")

DR.list <- plyr::llply(es.list, function(x) {
  x %>% as.data.frame %>% rownames_to_column()
  })
DR.list <- do.call(rbind, DR.list)
DR.Lambda.100Trees <- ggiraphExtra::summarySE(DR.list, measurevar = ".", groupvars = "rowname")
DR.Lambda.100Trees <- DR.Lambda.100Trees %>% rename(DR.Lambda.100Trees, DR = ., TipLabel = rowname)

all.trees.join <- left_join(DR.Lambda.100Trees %>% select(TipLabel, DR), BAMM.Lambda.100Trees %>% select(mean.lambda, TipLabel), by = "TipLabel")

BAMM.Lambda.100Trees %>% ggplot(aes(x = DR.Lambda.100Trees$DR, y = mean.lambda, fill = mean.lambda))+
  geom_point(shape = 21)+
  geom_errorbar(aes(ymin = mean.lambda - 1.96*sd, ymax = mean.lambda + 1.96*sd), size = 0.02)+
  geom_errorbarh(aes(xmin = DR.Lambda.100Trees$DR - 1.96*DR.Lambda.100Trees$sd, xmax = DR.Lambda.100Trees$DR + 1.96*DR.Lambda.100Trees$sd), size = 0.02)+
  theme_minimal()+
  ylab("100 Trees Speciation Rate (BAMM)")+
  xlab("100 Trees Speciation Rate (DR)")+
  scale_fill_viridis_c()

# mean(DR.Lambda.100Trees$sd)
# mean(BAMM.Lambda.100Trees$sd)

Figure S14: Speciation Rate means from 100 trees using either the DR statistic (x axis) or BAMM (y axis). While there is a clear correlation (r = 0.75) there is variability with BAMM showing less heterogeneity. 95 % CIs are plotted for both axis from the 100 trees. Each point represents a species (n = 5,965).

restricted.data %>% ggplot(aes(x = MCC.DR, y = MCC.BAMM.model.data$mean.lambda, fill = MCC.DR))+
  geom_point(shape = 21)+
  theme_minimal()+
  geom_smooth(method = "lm")+
  ylab("MCC Speciation Rate (BAMM)")+
  xlab("MCC Speciation Rate (DR)")+
  scale_fill_viridis_c()

Figure S15: Similar to Figure S14, there is a correlation (r = 0.68) between \(\lambda_{DR}\) and \(\lambda_{BAMM}\) with BAMM results showing less heterogeneity. Each point represents a species (n = 5,965).

R Session information

sessionInfo() %>% pander()

R version 3.3.1 (2016-06-21)

**Platform:** x86_64-apple-darwin13.4.0 (64-bit)

locale: en_AU.UTF-8||en_AU.UTF-8||en_AU.UTF-8||C||en_AU.UTF-8||en_AU.UTF-8

attached base packages: parallel, grid, stats, graphics, grDevices, utils, datasets, methods and base

other attached packages: bindrcpp(v.0.2.2), ggtree(v.1.6.11), EBImage(v.4.16.0), phylopath(v.1.0.1.9000), ggraph(v.1.0.1.9999), stargazer(v.5.2.2), tibble(v.1.3.4), RColorBrewer(v.1.1-2), kableExtra(v.0.7.0), phangorn(v.2.4.0), HDInterval(v.0.2.0), MuMIn(v.1.42.1), coda(v.0.19-1), car(v.3.0-2), carData(v.3.0-1), BAMMtools(v.2.1.6), dplyr(v.0.7.6), ggExtra(v.0.8), reshape2(v.1.4.3), purrr(v.0.2.4), caper(v.1.0.1), mvtnorm(v.1.0-6), MASS(v.7.3-50), ggridges(v.0.5.0), brms(v.2.4.0), Rcpp(v.0.12.18), phytools(v.0.6-44), maps(v.3.2.0), gridExtra(v.2.3), geiger(v.2.0.6), lme4(v.1.1-18-1), Matrix(v.1.2-14), repmis(v.0.5), diversitree(v.0.9-10), ape(v.5.1), mgcv(v.1.8-24), nlme(v.3.1-128), tidyr(v.0.7.2), knitr(v.1.20), pander(v.0.6.2) and ggplot2(v.3.0.0.9000)

loaded via a namespace (and not attached): R.utils(v.2.7.0), tidyselect(v.0.2.4), htmlwidgets(v.1.2), combinat(v.0.0-8), munsell(v.0.5.0), animation(v.2.5), codetools(v.0.2-14), units(v.0.5-1), DT(v.0.4), miniUI(v.0.1.1.1), ggiraphExtra(v.0.2.9.1), withr(v.2.1.2), Brobdingnag(v.1.2-6), colorspace(v.1.3-2), uuid(v.0.1-2), stats4(v.3.3.1), officer(v.0.2.1), bayesplot(v.1.6.0), labeling(v.0.3), rstan(v.2.17.3), mnormt(v.1.5-5), bridgesampling(v.0.5-2), rprojroot(v.1.3-2), clusterGeneration(v.1.3.4), R6(v.2.2.2), markdown(v.0.8), locfit(v.1.5-9.1), ggiraph(v.0.4.2), bitops(v.1.0-6), assertthat(v.0.2.0), promises(v.1.0.1), scales(v.1.0.0), gtable(v.0.2.0), tidygraph(v.1.1.0), rlang(v.0.2.2), scatterplot3d(v.0.3-41), splines(v.3.3.1), lazyeval(v.0.2.1), mycor(v.0.1), broom(v.0.5.0), inline(v.0.3.15), prediction(v.0.2.0), yaml(v.2.2.0), abind(v.1.4-5), threejs(v.0.3.1), crosstalk(v.1.0.0), backports(v.1.1.2), httpuv(v.1.4.5), rsconnect(v.0.8.8), tools(v.3.3.1), psych(v.1.8.4), gplots(v.3.0.1), BiocGenerics(v.0.20.0), plyr(v.1.8.4), base64enc(v.0.1-3), viridis(v.0.5.1), deSolve(v.1.20), zoo(v.1.8-3), haven(v.1.1.2), ggrepel(v.0.8.0), magrittr(v.1.5), data.table(v.1.11.4), openxlsx(v.4.1.0), colourpicker(v.1.0), sjmisc(v.2.6.3), R.cache(v.0.13.0), matrixStats(v.0.54.0), hms(v.0.4.2), shinyjs(v.1.0), mime(v.0.5), evaluate(v.0.10.1), fftwtools(v.0.9-8), xtable(v.1.8-2), shinystan(v.2.4.0), rio(v.0.5.10), jpeg(v.0.1-8), readxl(v.1.0.0), rstantools(v.1.5.1), crayon(v.1.3.4), KernSmooth(v.2.23-15), minqa(v.1.2.4), R.oo(v.1.22.0), StanHeaders(v.2.17.2), htmltools(v.0.3.6), later(v.0.7.3), tiff(v.0.1-5), udunits2(v.0.13), expm(v.0.999-2), sjlabelled(v.1.0.5), tweenr(v.0.1.5), ppcor(v.1.1), subplex(v.1.4-1), readr(v.1.1.1), cli(v.1.0.0), quadprog(v.1.5-5), R.methodsS3(v.1.7.1), gdata(v.2.18.0), bindr(v.0.1.1), igraph(v.1.1.2), forcats(v.0.3.0), pkgconfig(v.2.0.2), numDeriv(v.2016.8-1), foreign(v.0.8-66), xml2(v.1.2.0), dygraphs(v.1.1.1.6), stringdist(v.0.9.5.1), rvg(v.0.1.7), rvest(v.0.3.2), snakecase(v.0.7.1), stringr(v.1.2.0), digest(v.0.6.16), msm(v.1.6.6), rmarkdown(v.1.10), cellranger(v.1.1.0), fastmatch(v.1.1-0), gdtools(v.0.1.6), curl(v.3.2), shiny(v.1.1.0), gtools(v.3.8.1), nloptr(v.1.0.4), jsonlite(v.1.5), viridisLite(v.0.3.0), lattice(v.0.20-33), loo(v.2.0.0), httr(v.1.3.1), plotrix(v.3.7-3), survival(v.2.39-4), glue(v.1.3.0), xts(v.0.11-0), zip(v.1.0.0), png(v.0.1-7), shinythemes(v.1.1.1), ggforce(v.0.1.1), stringi(v.1.1.6) and caTools(v.1.17.1)

References

Armenta, J. K., P. O. Dunn, and L. A. Whittingham. 2008. Quantifying avian sexual dichromatism: A comparison of methods. Journal of Experimental Biology 211:2423. Journal Article.

BirdLife International and Handbook of the Birds of the World. 2017. Bird species distribution maps of the world. http://datazone.birdlife.org/species/requestdis.

Dale, J., C. J. Dey, K. Delhey, B. Kempenaers, and M. Valcu. 2015. The effects of life history and sexual selection on male and female plumage colouration. Nature 527:367–370. Journal Article.

Del Hoyo, J., A. Elliott, and D. Christie. 2011. Handbook of the birds of the world (Vols. 8-16). Book, Lynx Edicions 2003-2011.

Harvey, M. G., G. F. Seeholzer, B. T. Smith, D. L. Rabosky, A. M. Cuervo, and R. T. Brumfield. 2017. Positive association between population genetic differentiation and speciation rates in new world birds. Proceedings of the National Academy of Sciences 114:6328–6333.

Harvey Michael, G., L. Rabosky Daniel, and N. Cooper. 2017. Continuous traits and speciation rates: Alternatives to state‐dependent diversification models. Methods in Ecology and Evolution 9:984–993. Journal Article.

Jetz, W., G. H. Thomas, J. B. Joy, K. Hartmann, and A. O. Mooers. 2012. The global diversity of birds in space and time. Nature 491:444–448. Journal Article.

Quintero, I., and W. Jetz. 2018. Global elevational diversity and diversification of birds. Nature 555:246.

Schliep, K. P. 2010. Phangorn: Phylogenetic analysis in r. Bioinformatics 27:592–593.